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=  measurement  error  vector 
=  covariance  matrix  of  measurement  error  vector 
=  expected  value,  mean 
=  generalized  first  Cartesian  coordinate 
=  covariance  matrix  of  satellite  state 
oi  =  standard  deviation  of  ith  element  in  9  or  E 
X  =  covariance  matrix  of  estimated  parameters 
X  =  covariance  matrix  of  a  priori  parameter  estimates 
o  =  mass  fraction  of  all  mascons  to  lunar  mass 

(fe  =  latitude,  regular  or  geodetic,  of  site  or  mascon  i 

X  =  distance  ratio  between  the  lunar  maximum  radius  and  the 
lunar  satellite  radius 

¥,•  =  effect  of  central  body  attractions  of  body  i 
co  =  argument  of  perilune 
£1  =  right  ascension  of  the  ascending  node 
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Operators 

E{ }  =  expectation  operator 

e>q?U  =  exponential  operator 

ln[  ]  =  natural  logarithm  operator 
unit[  ]  =  unit  magnitude  operator  on  vector 

I  I  =  determinant  of  a  matrix,  magnitude  of  a  vector 
•  *  dot  product  operator 
X  =  cross  product  operator 

Subscripts 

b  =  lunar  satellite 
c  =  earth-moon  barycenter 
c.m.  =  center  of  mass 

CT  =  cross  track  direction 
DR  =  down  range  direction 
e  =  earth 

g  =  gravitational  parameters 
I  -  inertial  frame 
k  =  mascon 
1  =  moon 

L  =  local  vertical,  local  horizontal  frame 
m  =  order  (Legendre  functions  and  tesseral  coefficients) 
n  =  degree  (Legendre  functions  and  harmonic  coefficients) 
p  =  planet 

r  =  receive  time 

r. s.  =  receiving  site 

s  =  sun,  source  (signal),  or  send  time 

s. s.  =  sending  site 

t  =  tesseral  harmonics,  or  transpond/reflection  time 

ts.  =  transponding  or  reflection  site 
VT  =  vertical  direction 
z  =  zonal  harmonics 

0  »  initial  time 

i  =  ith  element  of  a  vector,  ith  column  vector  of  a  matrix 
ij  =  ijth  element  of  a  matrix 
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Superscripts 

R  =  position 

T  =  transpose  operator 

V  =  velocity 

-1  =  inverse 

=  vector 

=  vector  of  unit  length 
=  normalized  harmonic  coefficient 
=  first  derivative  with  respect  to  time 
'  =  first  derivative  with  respect  to  an  argument 
*  =  second  derivative  with  respect  to  an  argument 


Acronyms 

CT 

DSN 

GPS 

GSFC 

IAU 

JD 

JPL 

NASA 

PDI 

PEP 

PEP-D 

rms 

rss 

SAO 

UTC 


Coordinate  Time 
Deep  Space  Network 
Global  Positioning  System 
Goddard  Space  Flight  Center 
International  Astronomical  Union 
Julian  Date 

Jet  Propulsion  Laboratory 
National  Aeronautics  and  Space  Administration 
Powered  Descent  Initiation 
Planetary  Ephemeris  Program 
Draper  modified  Planetary  Ephemeris  Program 
root  mean  square  (square  root  of  the  mean  of  squared 
values) 

root  sum  square  (square  root  of  die  sum  of  squared  values) 
Smithsonian  Astrophysical  Observatory 

Universal  Coordinated  Time  (Universal  Temps  Coordin6e 
from  the  Bureau  International  de  l'Heure  in  France) 
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Chapter  One 


Introduction  and  Summary 


1.1  Background 

On  July  20, 1989,  die  20^  anniversary  of  the  first  Apollo  moon  landing. 
President  George  Bush  challenged  the  nation  to  undertake  an  ambitious 
course  of  human  space  exploration.  After  establishing  a  manned  presence  in 
earth  orbit  with  the  Space  Station  Freedom  in  the  1990’s,  the  President 
proposed  that  the  U.S.  return  to  the  moon,  and  return  to  stay.  From  this 
lunar  basing  point,  the  U.S.  could  continue  human  exploration  of  our  solar 
system  by  undertaking  a  manned  mission  to  Mars. 

The  establishment  of  a  lunar  base  will  result  in  significant  lunar  traffic 
to  build,  supply,  and  resupply  this  facility.  This  increased  traffic  will  require  a 
lunar  navigation  system.  As  the  nation  prepares  its  return  to  the  moon,  it 
will  have  to  decide  whether  this  navigation  system  should  be  earth-based, 
vehicle-based  or  lunar-based. 

Our  initial  voyages  to  the  moon  primarily  depended  upon  earth-based 
navigation  systems,  although  the  manned  missions  had  some  on-board 
capability.  An  earth-based  method  could  be  adopted  for  future  lunar  travel, 
but  NASA's  Deep  Space  Network  (DSN)  tracking  is  manpower  intensive, 
costly,  and  is  not  suitable  for  high  traffic  rates.  Additionally,  earth-based 
navigation  can  only  track  vehicles  on  the  lunar  near  side.  Adopting  an  earth- 
based  navigation  system  would  be  an  impractical  stepping  stone  for  human 
exploration  of  Mars  and  the  solar  system. 

A  vehicle-based  navigation  system  could  be  developed  to  support 
future  lunar  traffic.  Vehicle-based  navigation  has  become  practical  because  of 
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advances  in  inertial  navigation  equipment  and  on-board  computing 
capabilities.  Inertial  navigation  would  be  limited  by  our  knowledge  of  the  ds- 
lunar  environment,  principally  our  knowledge  of  the  moon's  gravity  field. 

The  limitations  of  earth-based  navigation  and  the  high  accuracy 
requirements  of  certain  mission  phases  (principally  landing)  may  require  the 
development  of  a  lunar-based  navigation  system.  A  simple  lunar  navigation 
system  similar  to  the  earth's  Global  Positioning  System  (GPS)  would  handle 
high  traffic  rates  and  would  provide  accurate  lunar  far-side  navigation.  If 
high  lunar  traffic  rates  are  achieved,  then  a  lunar-based  system  could  provide 
a  higher  accuracy  system  alternative  to  vehicle-based  navigation. 

Many  low-altitude  mission  phases  will  require  accurate  knowledge  of 
the  moon's  gravitational  field,  especially  its  far-side  characteristics.  Spherical 
harmonic  models  are  typically  used  to  model  the  gravitational  field  near 
celestial  bodies.  Finite  expansion  spherical  harmonic  models,  however,  do 
not  accurately  model  the  low  altitude  gravity  field  of  the  moon.  This  is 
because  the  moon's  gravitational  field  contains  significant  anomalies, 
discovered  in  the  late  1960's  by  scientists  at  NASA's  Jet  Propulsion  Laboratory 
(JPL)  [35],  making  such  models  inefficient.  From  earth-based  lunar  tracking 
data,  the  scientists  developed  a  lunar  gravity  model.  This  model  was  then 
compared  to  topographical  images  of  the  moon  and  revealed  mass 
concentrations  around  the  ringed  maria.  These  mass  concentrations  or 
"mascons"  exhibit  very  high  frequency  gravitational  behavior  and  therefore 
require  a  very  high  number  of  terms  in  the  spherical  harmonic  expansion  to 
model  this  behavior.  Since  their  discovery,  scientists  have  postulated 
different  models  to  account  for  the  lunar  mascon  phenomenon,  since 
expanding  the  spherical  harmonic  model  to  high  degree  and  order  was 
computationally  impractical.  Chapter  Two  surveys  the  spherical  harmonic 
and  several  other  gravitational  field  modeling  techniques  in  more  detail. 


1.2  Motivation 

Since  the  real  lunar  gravitational  field  is  difficult  to  accurately  model, 
this  thesis  will  study  the  implications  of  modeling  errors.  An  inaccurate 
model  of  the  lunar  gravity  field  will  result  in  the  growth  of  navigation  errors. 
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Mismodeled  acceleration  forces  result  in  both  velocity  and  position  errors. 
Since  gravitational  acceleration  is  a  function  of  position  position  errors  will 
lead  to  increased  acceleration  errors,  further  increasing  the  velocity  and 
position  errors.  This  error  propagation  may  or  may  not  be  critical  depending 
upon  the  magnitude  of  the  errors,  the  navigation  system's  ability  to  measure 
them,  and  the  mission  phase  accuracy  requirements. 

Specifically,  an  inaccurate  gravity  model  will  significantly  affect  any 
landing  maneuvers  with  strict  accuracy  requirements.  Unmanned  cargo 
missions  to  resupply  a  lunar  base  will  be  particularly  vulnerable  to  errors 
from  a  mismodeled  gravity  field.  Since  there  is  no  appreciable  lunar 
atmosphere,  gravity  forces  dominate  a  vehicle's  descent  to  the  moon's 
surface.  Since  the  force  of  gravity  is  inversely  proportional  to  the  square  of 
distance,  navigation  errors  due  to  a  mismodeled  gravity  field  increase  as  the 
vehicle  descends  to  the  surface.  Avoiding  unacceptable  landing  errors  will 
depend  upon  an  accurate  determination  of  the  lunar  gravity  field. 

The  scientific  community  is  also  interested  in  developing  a  more 
precise  model  of  the  lunar  gravity  field.  A  better  model  can  improve 
knowledge  of  the  moon's  composition  and  internal  structure.  Models  of 
different  elements,  their  densities,  and  their  distribution  within  the  moon's 
interior  could  be  developed  to  match  the  observed  gravitational  field. 
Gravitational  models  may  also  help  to  determine  the  selenological  thermal 
and  tectonic  history.  The  discovery  of  lunar  mas  cons  has  also  led  to  scientific 
speculation  about  how  mass  concentrations  formed  in  these  shallow  seas. 
The  scientific  community  hopes  that  a  better  understanding  of  the 
gravitational  field  around  the  ringed  maria  and  other  lunar  surface  features 
will  help  to  determine  the  origin  of  these  features  [2, 35]. 

The  purpose  of  this  thesis  is  to  determine  the  feasibility  of  using  a 
spherical  harmonic  lunar  gravitational  model,  based  on  observations  of  a 
near-circular  polar  satellite,  to  predict  low  altitude  lunar  orbits  globally. 
Rather  than  attempting  to  develop  a  more  precise  lunar  gravitational  field 
model,  this  thesis  investigates  measurement  types  and  satellite  orbits  that  can 
be  vised  to  develop  gravity  field  models.  Each  measurement  type  will  have 
advantages  and  disadvantages  in  terms  of  cost,  schedule,  and  accuracy.  This 
thesis  investigates  each  different  method's  ability  to  estimate  a  lunar 
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gravitational  potential  model  and  the  model's  accuracy  in  predicting  future 
lunar  orbits.  By  analyzing  the  capabilities  of  different  sensing  methods,  this 
investigation  will  allow  NASA  to  plan  unmanned  lunar  precursor  missions 
to  extract  die  best  lunar  gravity  field  information. 

13  Initial  Lunar  Gravitational  Sensing  Method 

Current  lunar  gravitation  field  models  are  based  upon  earth-based 
tracking  data  from  the  Lunar  Orbiter  program  of  the  1960's  (Figure  1.3-1). 
These  unmanned  Apollo  precursor  missions  provided  photographic  imaging 
and  gravitational  mapping  of  the  moon.  Apollo  lunar  navigation  was  based 
upon  Lunar  Orbiter's  gravity  field  mapping.  In  addition  to  the  Lunar  Orbiter 
missions,  tracking  data  from  Apollo  missions  and  some  Soviet  lunar 
missions  are  included  in  current  gravity  field  models  [3, 12, 19, 33, 41, 47] 

Lunar  Orbiter  gravitational  mapping  missions  utilized  Doppler 
measurements  of  radio  tracking  signals.  Lunar  Orbiter  spacecraft  were  tracked 
by  NASA's  Deep  Space  Network  (DSN)  across  the  near  side  of  the  moon  [36]. 
A  DSN  tracking  station  sent  a  continuous  wave  S-Band  frequency  to  the 
spacecraft.  The  spacecraft  received  this  Doppler-shifted  signal  and  re¬ 
transmitted  it  to  earth.  The  tracking  station  received  this  signal,  Doppler- 
shifted  once  again  in  frequency.  The  tracking  station  used  this  signal  to 
calculate  the  relative  velocity  between  the  spacecraft  and  tracking  station.  The 
relative  velocity  observed  was  then  combined  with  position  tracking  data  to 
estimate  the  lunar  gravity  field  using  methods  similar  to  those  discussed  in 
Chapter  Five. 


Figure  13-1:  Lunar  Orbiter  Earth-Moon  Geometry 


Current  lunar  gravitation  field  model  accuracy  is  limited  by  the 
amount  of  lunar  orbital  tracking  data  available.  The  Lunar  Orbiter  and 
Apollo  missions  were  mostly  low  inclination  missions.  Since  most  missions 
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flew  about  the  lunar  equator,  derived  gravity  field  models  emphasize  the 
effects  of  anomalies  near  the  equator.  Only  a  fraction  of  the  moon's  surface 
was  covered  by  these  missions  [32],  and  therefore  gravity  field  data  is  lacking 
for  lunar  polar  regions.  Additionally,  this  earth-based  gravity  mapping 
method  was  limited  to  observations  of  near-side  lunar  spacecraft  passes. 
Gravitational  disturbances  on  the  far  side  were  only  determined  by  their 
integrated  effects  on  satellite  position  and  velocity  from  the  end  of  one  near¬ 
side  pass  to  the  beginning  of  the  next  pass.  Thus  current  lunar  gravitational 
Held  models  do  not  provide  very  meaningful  information  about  the  lunar  far 
side. 


1.4  Proposed  Lunar  Gravitational  Sensing  Methods 

Any  future  lunar  gravitational  field  sensing  system  will  have  to  greatly 
improve  our  knowledge  of  the  moon's  gravitational  field  to  justify  the 
mission's  cost.  To  achieve  this  improvement  in  accuracy,  the  system  will 
have  to  address  the  current  model's  limitations.  The  motion  of  an  orbiting 
body  should  be  sensed  without  any  orbital  maneuvering  which  disturbs  the 
estimation  solution.1  Thus  it  is  desirable  to  select  orbits  which  are  stable  for  at 
least  one  lunar  orbit  to  avoid  re-boost  maneuvering.  A  high  inclination, 
preferably  polar,  lunar  orbiter  would  allow  observations  of  satellite 
accelerations  over  the  moon's  entire  surface  as  the  moon  rotates  under  the 
orbital  plane.  A  dual  orbiter  sensing  scheme  would  be  better  because  it  would 
allow  lunar  far-side  accelerations  to  be  observed.  Better  still  would  be  a 
sensing  scheme  observing  the  motion  of  several  satellites  in  different 
inclinations  than  those  available  in  Apollo-era  lunar  missions. 

NASA  is  considering  two  different  sensing  schemes.  NASA's  Jet 
Propulsion  Laboratory  has  proposed  a  dual  orbiter  scheme  which  uses  radio- 
based  Doppler  observations  to  sense  the  moon's  gravitational  field  effects  [40]. 
NASA's  Goddard  Space  Flight  Center  has  proposed  a  co-orbital  scheme  which 


1  Gravitational  accelerations  experienced  by  an  orbiter  are  not  measured  directly.  Methods  to 
measure  the  accelerations  due  to  gravity  therefore  use  external  observations  of  a  body's  motion. 
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uses  a  laser-based  system  which  makes  both  ranging  and  Doppler 
observations  [2]. 


The  dual  orbiter  sensing  scheme  uses  a  low  altitude,  circular  polar 
satellite  and  a  high  apolune  polar  elliptical  satellite  (Figure  1.4-1).  The 
elliptical  satellite  orbit  is  positioned  so  that  apolune  initially  occurs  on  the  far 
side  of  the  moon.  This  increases  the  duration  of  lunar  far-side  viewing.  The 
orbit  is  also  skewed  such  that  apolune  is  outside  of  the  earth  occultation  zone 
which  increases  the  time  that  the  elliptical  satellite  is  within  the  line  of  sight 
of  earth  tracking  stations.  From  these  orbits  the  relative  velocity  between  the 
two  spacecraft  can  be  measured  using  either  the  bent  pipe  or  the  satellite 
bounce  methods  described  below. 


Low-Altitude 
Spacecraft 
Orbit 


Two-Way  Coherent  Doppler 


Elliptical  Satellite 
Orbit 


Earth 

Occultation 

Zone 


Two-Way 

Coherent 

Doppler, 


Figure  1.4-1:  Dual  Orbiter  Sensing  Method 


The  bent  pipe  method  uses  a  four-way  coherent  Doppler  scheme  in 
which  a  high  frequency  is  generated  by  an  atomic  clock  at  a  DSN  tracking  site 
and  transmitted  to  the  elliptical  "viewing"  satellite.  The  "viewing"  satellite 
uses  the  Doppler-shifted  received  signal  to  generate  a  lower  frequency  signal 
which  it  transmits  to  the  circular  "gravity  sensing"  spacecraft.  This  spacecraft 
receives  the  Doppler-shifted  signal  and  re-transmits  it  to  the  "viewing" 
satellite.  The  "viewing"  satellite  modulates  the  received  Doppler-shifted 
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frequency  signal  from  the  "sensing7'  spacecraft  onto  the  frequency  signal 
received  from  the  tracking  station  and  transmits  this  signal  back  to  earth.  The 
signal  received  by  the  tracking  station  is  then  processed  to  retrieve  both  the 
relative  velocities  between  the  tracking  station  and  "viewing"  satellite  and 
between  the  "viewing"  and  "gravity  sensing"  satellites. 


The  satellite  bounce  method  uses  a  two-way  coherent  Doppler  scheme 
between  the  two  spacecraft.  The  circular  "gravity  sensing"  satellite  generates  a 
continuous  wave  frequency  signal  for  transmission  to  the  elliptical  "viewing" 
satellite.  This  satellite  shifts  the  signal's  frequency  for  transmission  back  to 
the  first  one.  The  receiving  spacecraft  extracts  the  Doppler  shift  from  the 
signal,  records  it  and  transmits  it  to  earth  when  in  view  of  an  earth  tracking 
station.  Coherent  Doppler  links  between  earth  tracking  stations  and  either 
satellite  are  used  to  aid  in  the  estimation  of  the  lunar  gravitational  field. 


Earth-based 
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Figure  1.4-2:  Co-orbital  Sensing  Method 


NASA's  Goddard  Space  Flight  Center  has  proposed  a  co-orbital  scheme 
in  which  a  satellite  in  a  circular  polar  orbit  ejects  a  subsatellite  in  the  same 
orbit  (Figure  1.4-2).  Both  spacecraft  are  affected  by  lunar  gravitational 
perturbations,  so  both  are  "sensing"  vehicles  and  a  laser  system  measures 
their  relative  motion.  The  satellite  contains  the  sensing  equipment  and  the 
subsatellite  is  a  passive  reflector  for  the  satellite's  emitted  laser  beams.  The 
satellite's  laser  transmits  a  light  beam  toward  the  subsatellite.  This  beam 
reflects  off  of  the  subsatellite  back  to  the  satellite.  A  high  accuracy  ranging 
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measurement  between  the  satellite  and  subsatellite  is  made  by  determining 
the  travel  time  of  a  transmitted  sub-carrier  pulse  signal.  Lesser  accuracy 
Doppler  relative  velocity  measurements  are  made  from  the  frequency 
shifting  of  the  transmitted  laser  signal.  Observation  data  is  stored  and 
transmitted  to  earth  at  regular  intervals. 


1.5  Simulation  Tools 

The  primary  tool  used  to  accomplish  the  goals  of  this  thesis  is  the 
Planetary  Ephemeris  Program  (PEP),  a  FORTRAN  computer  program 
obtained  from  the  Smithsonian  Astrophysical  Observatory  (SAO)  and 
executed  on  Sim  workstations  at  the  Charles  Stark  Draper  Laboratory.  The 
Smithsonian  version  of  PEP  has  most  of  the  capabilities  needed  for  the 
analyses  of  this  thesis.  Modifications,  coded  at  Draper  Laboratory,  have 
augmented  its  capabilities  for  this  thesis  research.  In  addition,  auxiliary 
software  has  been  developed  to  analyze  lunar  orbits,  die  observation  schemes, 
navigation  uncertainties,  and  estimated  lunar  gravitational  fields. 

Given  the  description .  of  a  body's  gravitational  field,  PEP  can 
numerically  integrate  a  satellite's  motion  about  that  body.  For  this  thesis,  PEP 
was  modified  to  accommodate  a  point  mass  (mascon)  gravity  model  in 
addition  to  the  spherical  harmonic  model.  The  techniques  used  in  PEP  for 
numerically  integrating  the  differential  equations  of  motion  for  a  lunar 
satellite  cure  described  in  Chapter  Three.  This  chapter  also  describes  the 
methods  PEP  uses  to  calculate  the  partial  derivatives  of  the  motion  with 
respect  to  orbital  initial  conditions,  gravity  harmonic  coefficients,  and  other 
parameters.  The  mascon  gravity  model  modifications  coded  in  PEP  are  also 
covered  in  Chapter  Three. 

The  Planetary  Ephemeris  Program  (PEP)  can  also  process  the 
astronomical  observations  generated  by  the  various  lunar  gravitational 
sensing  methods.  The  many  different  types  of  observations  that  can  be 
processed  in  PEP  are  described  in  Chapter  Four.  This  chapter  also  discusses 
how  PEP  generates  simulated  observations  with  a  truth  model  (mascon  and 
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spherical  harmonic  expansion)  and  then  estimates  theoretical  values  of  the 
observations  for  another  model  (spherical  harmonic  expansion)  to  fit  those 
"truth  model"  observations.  Least  squares  maximum  likelihood  estimation 
and  prediction  uncertainty  propagation  techniques  are  described  in  Chapter 
Five.  For  this  thesis/  noise  has  not  been  included  in  the  satellite  dynamics. 
Kalman  filter  and  system  identification  techniques  may  be  needed  when 
processing  real  observations,  because  of  the  noise  due  to  radiation  pressure, 
gas  leakage,  and  other  unmodeled  forces.  These  estimation  techniques  are 
also  described  in  Chapter  Five. 


1.6  Methodology 

The  focus  of  this  thesis  is  the  estimation  of  a  lunar  gravity  field  model 
based  on  various  measurement  types  and/or  orbital  geometries.  Each  unique 
measurement  type  and  orbital  geometry  combination  will  be  refered  to  as  a 
sensing  scheme.  The  standard  earth-based  orbiter  state  sensing  scheme  and 
the  proposed  dual  orbiter  bent  pipe  scheme  are  analyzed  in-depth. 
Additionally,  the  co-orbital  laser  ranging  scheme,  a  non-coplanar  bent  pipe 
scheme,  and  an  earth-based  interferometric  observation  scheme  are 
investigated  to  determine  whether  any  of  these  schemes  can  reduce  the 
parameter  correlation's  observed  during  gravitational  parameter  estimation. 

Since  the  true  lunar  gravitational  field  is  not  precisely  known,  a 
"truth"  model  was  developed  and  used  for  this  investigation,  as  discussed  in 
Section  6.4.  This  "truth"  model  combines  Bills  and  Ferrari’s  16  X  16  lunar 
harmonic  model  [12]  up  to  degree  and  order  five,  along  with  78  point  masses 
distributed  below  the  lunar  surface  to  simulate  the  behavior  of  mass 
concentrations.  This  truth  model  was  used  to  simulate  observations  for  the 
various  sensing  schemes. 

For  each  sensing  scheme,  the  coefficients  in  spherical  harmonic  fit 
models  of  degree  and  order  8  and  12  were  estimated  to  optimally  represent 
the  "true"  gravity  field  model  by  fitting  to  the  "truth"  model  observations. 
Using  first  guesses  for  the  fit  model  coefficients  and  satellite  initial  osculating 
orbital  elements,  the  equations  of  motion  and  the  equations  for  the  partial 
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derivatives  of  the  motion  with  respect  to  these  quantities  were  numerically 
integrated.  From  these  numerical  integrations,  theoretical  values  of  the 
observations  were  determined.  The  observation  residuals  (difference 
between  the  "truth"  model  observations  and  fit  model's  theoretical 
observation  values)  and  the  observation  partial  derivatives  were  computed 
and  used  to  obtain  parameter  adjustments  to  the  first  guesses  for  the  fit  model 
coefficients  and  satellite  initial  osculating  orbital  elements.  This  process  was 
then  repeated  until  either  the  method  converged  upon  a  solution  or  no 
solution  could  be  determined.  Chapter  Six  discusses  the  implementation  of 
the  estimation  process  for  various  sensing  methods  and  certain  test  cases. 

The  estimated  lunar  gravitational  field  models  were  analyzed  to 
determine  their  accuracy  relative  to  the  "truth".  Navigation  errors 
propagated  in  one  lunar  revolution  were  used  to  determine  the  estimated 
gravitational  field  model’s  accuracy.  Two  types  of  lunar  orbits  were  analyzed: 
a  15°  inclination,  100  km  altitude,  near-circular  orbit  and  a  lunar  landing 
from  a  5°  inclination,  200  km  near-circular  parking  orbit.  These  orbits  were 
propagated  using  both  the  "truth"  and  estimated  models.  The  analysis  of  the 
estimated  gravitational  field  models  is  discussed  in  detail  in  Chapter  Seven. 
In  addition,  this  chapter  discusses  the  attempts  made  to  break  the  high 
parameter  correlations  which  were  discovered  during  the  estimation  process. 

For  the  near-circular  orbit,  the  position  and  velocity  errors  between  the 
"truth"  and  estimated  gravity  fields  were  used  to  quantify  the  estimated 
model’s  accuracy  and  thus  the  sensing  method's  capability.  The  gravitational 
parameter  covariance  matrix,  determined  during  the  estimation  process,  was 
used  with  the  estimated  field's  orbit  propagation  with  partial  derivatives  to 
predict  position  and  velocity  uncertainties.  These  predicted  uncertainties 
were  then  compared  to  the  true  state  errors  between  the  two  models'  orbits. 
In  real  gravity  field  missions,  the  true  gravity  field  will  not  be  available  for 
comparison  with  the  estimate,  so  it  is  useful  to  understand  the  relation 
between  these  two  analyses. 

For  the  lunar  landing  maneuver,  the  estimated  gravity  field  model  was 
used  to  determine  the  deorbit  bum  and  the  selenographic  position  for 
Powered  Descent  Initiation  (PDI).  The  spacecraft's  circular  parking  orbit  was 
numerically  integrated  until  the  appropriate  time  for  the  deorbit  bum.  The 
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spacecraft's  velocity  at  bum  time  was  then  adjusted  to  simulate  the  deorbit 
maneuver  and  the  spacecraft's  elliptical  transfer  orbit  was  then  integrated. 
The  spacecraft's  position  upon  reaching  the  powered  descent  stage  of  the 
mission  was  then  determined.  The  "true"  PDI  point  was  then  compared  to 
the  estimated  model  target  PDI  point  to  determine  the  model's  accuracy  for 
planning  lunar  landing  maneuvers. 


1.7  Summary  of  Results 

The  analyses  of  the  8  X  8  estimated  fit  models  clearly  demonstrated  that 
lunar  far-side  observations  are  required  in  the  accurate  estimation  of  the 
lunar  gravity  field.  For  both  the  lunar  landing  maneuver  and  the  satellite 
state  uncertainty  prediction,  the  observation  technique  which  included  lunar 
far-side  observations  in  addition  to  earth-based  near-side  observations 
produced  a  much  more  accurate  lunar  gravity  fit  model.  This  estimated 
model  planned  a  lunar  deorbit  maneuver  4.3  times  more  accurately  than  the 
model  based  on  earth-based  observations  alone.  The  earth-based  observation 
fit  model  also  predicted  state  uncertainties  four  times  larger  than  its 
counterpart.  The  earth-  and  satellite-based  fit  model  produced  state  errors  for 
the  single  orbit  that  were  again  roughly  a  quarter  of  the  earth-based  fit 
model's  errors. 

The  lunar  navigation  analyses  also  demonstrate  that  the  eighth  degree 
and  order  spherical  harmonic  expansion  fit  models  were  unable  to  adequately 
model  ttie  lunar  mascons  included  in  the  lunar  gravitational  "truth"  model. 
In  the  best  cat  e,  the  8  x  8  fit  model  predicted  single  orbit  ahead  uncertainties  of 
close  to  three  quarters  of  a  kilometer  in  position  and  one  meter  per  second  in 
velocity.  The  orbit's  actual  state  errors  were  closer  to  three  kilometers  in 
position  and  two  and  a  half  meters  per  second  in  velocity.  Additionally,  the 
best  fit  model  produced  a  fifty-six  kilometer  position  error  for  the  lunar 
deorbit  maneuver.  These  results  are  discussed  in  further  detail  in  Chapter 
Eight,  which  also  recommends  several  subjects  related  to  this  thesis  which 
deserve  further  study. 
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Chapter  Two 


Gravity  Field  Models 


2.1  LaPlace's  Equation  and  die  External  Gravitational  Field 

According  to  Newton's  law  of  gravitation  two  particles  attract  each 
other  with  a  force,  acting  along  the  line  joining  them,  which  is  proportional 
to  the  product  of  their  masses  and  inversely  proportional  to  the  square  of  the 
distance  between  them  [10,  p.  95].  From  this  law,  the  gravitational  force  of  a 
body  of  mass  m2  acting  upon  a  body  of  mass  mi  can  be  mathematically 
represented  by  the  formula: 

f=  (2.1-1) 

FWiI 

where  G  is  the  universal  gravitational  constant,  and  r,  and  f2  are  the  position 
vectors  of  bodies  one  and  two  respectively. 

Unfortunately,  this  formula  is  only  appropriate  if  the  two  bodies  are 
point  masses,  or  behave  as  them.  Such  is  the  case  for  spherical  bodies  if 
density  is  a  function  of  radius  from  the  center  only.  The  point  mass  model 
also  provides  an  accurate  representation  of  the  gravitational  attraction  for 
widely  separated  bodies.  In  the  limit  of  large  distances,  gravitational  bodies 
tend  to  look  like  point  masses  so  that  mass  distribution  becomes 
unimportant. 

For  many  practical  applications,  the  attracting  body  cannot  be  modeled 
as  a  point  mass  and  a  different  mathematical  model  must  be  developed.  For 
the  case  in  which  the  attracted  body  is  small  compared  to  the  attracting  body 
and  the  attracting  body  is  an  arbitrary  distribution  of  mass  with  a  finite 
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dimension,  the  force  on  a  mass  m  located  at  position  vector  F  in  Figure  2.1-1 
produced  by  an  element  of  mass  dM  with  position  vector  R  is 


dF(f)  =  -GmdM 


(2.1-2) 


Figure  2.1-1:  Distributed  Mass  acting  upon  a  Point  Mass  m 

Integration  over  the  entire  volume  of  the  distributed  body  will  produce  the 
gravitational  force  on  the  mass  m  of  the  attracting  body  of  total  mass  M.  This 
force  can  then  be  represented  by  a  scalar  potential  U,  such  that  the 
gravitational  force  on  a  body  located  outside  of  the  attracting  body  may  be 
obtained  as  the  gradient  of  the  scalar  potential,  or 

F(r)  =-m  VU(r).  (2.1-3) 

There  is  a  sign  convention  discrepancy  between  some  of  the  references  used 
and  PEP  documentation  [8].  According  to  Kaula  [24],  physicists  define  the 
gradient  of  the  potential  field  as  in  Equation  (2.1-3),  whereas  astronomers 
define  the  same  gradient  with  a  change  in  sign.  The  formulas  in  this  thesis 
follow  the  former  convention  to  agree  with  PEP  documentation  and  software 
coding  modifications  explained  herein. 


The  scalar  gravitational  potential  U  can  be  written  in  terms  of  the  rectangular 
coordinate  system  of  Figure  2.1-1  and  the  mass  density  p  as  [14] 
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The  Lapladan  of  the  scalar  potential  U  in  rectangular  coordinates  is 

Taking  the  required  partial  derivatives  of  the  potential  U  from  (2.1-4)  yields 

-3 


V2U  =  -GM 


[(z-|)2+(y-7,)2+(z-C)2f  + 

[(*-{)* +(y-|?)*+(g-01| 
[(*-l)2+(y-'j)2+(z-?)2] 


x 


Since  the  bracketed  terms  cancel,  this  reduces  to 

V2UsO 


(2.1-6) 


(2.1-7) 


This  relationship  is  known  as  Laplace's  equation  and  applies  at  all  points 
outside  of  the  distributed  attracting  mass.  Its  solutions  are  called  harmonic 
functions.  Any  scalar  function,  U,  which  satisfies  Laplace's  equation  and  the 
far-field  boundary  condition  that  the  potential  approaches  0  as  1/r  can  be  used 
to  desdbe  the  gravitational  field  about  some  distributed  mass.  If  U  is  defined 
with  suffident  flexibility,  i.e.  an  infinite  number  of  orthogonal  terms  with 
undefined  constant  coeffidents,  then  U  can  be  tailored  to  describe  the 
gravitational  field  outside  of  any  arbitrarily  distributed  mass.  The  above 
method  of  deriving  (2.1.7)  is  based  upon  the  method  used  by  Kaula  [24],  Battin 
[10],  and  Comfort  [14]. 


2.2  Spherical  Harmonic  Expansion  for  the  Gravitational  Potential 

The  most  common  gravitational  potential  model  is  the  spherical 
harmonic  expansion.  This  expansion  can  be  derived  by  solving  Laplace's 
equation  in  spherical  harmonic  coordinates.  First,  the  rectangular  coordinate 
system  of  Figure  2.1-1  is  converted  to  spherical  coordinates  through  the 
transformation  below,  which  is  depicted  in  Figure  2.2-1.  If  the  center  of  the 
distributed  mass  is  selected  as  the  origin  of  this  coordinate  system,  the 
expansion  is  simplified. 
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x-n 
y  =  r  sin  6  cos  <p 
2  =  rsin0 


0  <  r  <  oo 
O<0<2tt 


Figure  2.2-1:  Spherical  Coordinate  System 
In  spherical  coordinates,  Laplace's  equation  becomes  [24] 


1 


cos0  d<f> 


3  (  .dV'' 


1  32  U 


(25-1) 


COS2  0  do2 


=  0  (2.2-2) 


One  common  method  of  deriving  the  solution  to  U  in  spherical 
coordinates  is  through  the  method  of  separation  of  variables,  maintaining  the 
boundary  condition  on  r.  This  method  can  be  found  in  Kaula  [24]  or  Comfort 
[14]  and  leads  to  the  solution 

u M,0)  =  Cosm0+S„„sinm0]  (2.2-3) 

»=0  m=0 

where  Cnm  and  Snm  are  constant  coefficients  and  Pnm  are  the  generalized 
Legendre  functions  of  degree  n  and  order  m.  Equation  (2.2-3)  is  the  complete 
real  solution  of  Laplace's  equation  in  spherical  coordinates.  There  are  other 
solutions  which  are  physically  inadmissible  since  they  cause  the  potential  to 
become  infinite  at  <J>  =  n/2,  3n/2.  These  solutions  involve  the  associated 
Legendre  functions  of  the  second  kind.  Another  method  of  deriving  the 
solution  to  U  is  to  expand  the  denominator  of  (2.1-4)  in  Legendre  functions. 
This  is  the  procedure  carried  out  by  Battin  [10]  and  Croopnick  [18]. 

Equation  (2.2-3)  is  the  generalized  solution  of  Laplace's  equation  and 
must  be  modified  for  the  gravitational  case.  Laplace's  equation  holds  true 
everywhere  outside  the  surface  of  the  gravitational  body.  The  harmonic 
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expansion,  however,  is  convergent  only  from  pi  to  «  in  the  radial  direction, 
where  pi  is  the  body's  (lunar)  maximum  radius,  usually  the  equatorial  radius. 
Furthermore,  pi  can  be  used  to  non-dimensionalize  the  Cnm  and  Snm 
coefficients  as  shown  in  Equation  (2.2-4).  Secondly,  U  is  multiplied  by  the 
gravitational  constant  and  the  total  mass  of  the  attracting  body  to  satisfy 
Equation  (2.1-3).  The  n=0,  m=0  term  is  the  point  mass  approximation  for  the 
attracting  body  and  can  be  separated  from  the  summation  for  ease  of 
computation.  The  Snm  coefficient  is  meaningless  for  m=0  because  sin(mO)  is 
always  zero.  Separating  the  m=0  terms  simplifies  the  summation.  The 
Pno(sin<)>)  Legendre  functions  then  become  the  standard  Pn(sin$)  Legendre 
polynomials.  When  the  m=0  terms  are  separated,  the  Cno  terms  are 
commonly  replaced  by  Jn  terms,  where  Jn  =  -Cno-  After  the  previous 
modifications,  the  spherical  harmonic  gravitational  potential,  U,  becomes 


U(r,0,fl)  = 


-GM 


-m 


Pn  (sin  <f>)  + 


X(t]  c°sm0+s„„ sinmfl] 

n=l  ^  T  '  m=l 


(2.2-4) 

Additionally,  if  the  attracting  body's  center  of  mass  is  the  origin  of  the 
coordinate  system,  then  the  first  degree  terms  are  identically  zero  [22], 

Ji  =0,  C„  =  0,  Sn=0  (2.2-5) 


and  the  summation  limits  can  begin  with  the  second  degree  terms. 


When  dealing  with  large  degree  models,  it  is  common  to  normalize 
the  Legendre  functions,  and  therefore  adjust  the  tesseral  coefficients  as  well. 
The  Legendre  functions  are  normalized  to  satisfy  the  equation 


j(P„m(sm<(>)cosmd)2ddd<l>  =  4it 
J  (P„m(sin  <p)smm6)2  dQd<p  =  An 


O<0<2  n 


{-%<(!><% 


(2.2-6) 


This  leads  to  the  following  relationship  between  the  unnormalized  and 
normalized  Legendre  functions  and  the  inverse  relationship  between  the 
unnormalized  and  normalized  coefficients: 
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(m  =  0) 


(2  2-7) 


The  Planetary  Ephemeris  Program  (PEP)  was  modified  to  allow 
normalized  Legendre  functions  to  be  used  for  higher  degree  harmonic 
models.  Normalized  coefficients  are  better  conditioned  for  PEP'S  floating 
point  parameter  estimation  algorithms,  especially  for  higher  degree  models. 
Appendix  A  discusses  the  Legendre  polynomials,  the  generalized  Legendre 
functions  and  the  recursive  formulas  developed  for  their  implementation. 
Although  use  of  the  normalized  Legendre  functions  was  considered  necessary 
for  this  thesis,  this  scaling  was  not  considered  necessary  for  the  Legendre 
polynomials.  Because  of  this  inconsistency  in  scaling,  the  equations  coded  in 
the  software  use  normalized  tesseral  coefficients  and  unnormalized  zonal 
coefficients. 


After  accounting  for  the  above-mentioned  modifications,  the  spherical 
harmonic  potential  equation  (2.2-4)  becomes 


U(r,M)  = 


-GM 


p«(sin^) 


S(t)  cos'n®  +  sm» smmfl] P„„(sin0) 

V  '  J  m_i 


(2.2-9) 


The  Planetary  Ephemeris  Program  uses  this  equation  with  finite  summation 
upper  limits.  The  coding  in  PEP  separates  the  central  body  term  from  the 
harmonic  summations. 
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2.2.1  Relation  between  Spherical  Harmonics  and  Moments  of  Inertia 


The  moments  of  inertia  about  a  solid  body  are  closely  related  to  the 
gravitational  potential  for  a  distributed  mass.  An  expansion  of  Equation  (2.1- 
4)  for  the  gravitational  potential  contains  the  terms  for  the  moments  of 
inertia. 

body  body 

lm=\j\(S2+Z1)dm  =  JjJ$?  *"  (2.21-D 

body  body 

I«=j]j(*2  +  '»2)*" 

body  body 


Equations  (2.2.1-1)  are  valid  for  any  arbitrary  rectangular  coordinate 
system  (£.11,0  originating  at  the  body's  center  of  mass.  The  spherical 
harmonic  potential  function  U  was  developed  from  the  (r,$,6)  coordinate 
system  of  Figure  2.2-1  which  corresponds  to  the  (x,y,z)  rectangular  coordinate 
system  of  Figure  2.1-1.  Converting  the  moments  of  inertia  to  this  (x,y,z) 
system,  the  second  degree  coefficients  in  the  gravitational  potential  function 
satisfy  the  following  relations  [24, 30, 32, 44J. 


_  K(i„+iw)-i= 

20 - riM 

r  la  c  _ 
pfM  p}M 


-  _  ” 
““  4p,2M 

S_  *v 
22”  7 

2 pfM 


(2.2.1-2) 


Since  the  lunar  moments  of  inertia  can  be  obtained  by  observing  the 
physical  librations  of  the  moon,  the  second  degree  coefficients  can  be 
determined  without  sending  spacecraft  to  the  moon.  If  the  (x,y,z)  axes  are 
principle  axes,  the  products  of  inertia  will  be  zero,  resulting  in 


C21  =  0,  S21  =  0,  S22  =  0 


(2.2.1-3) 


In  the  lunar  case,  the  principle  axes  coincide  with  the  x  axis  pointing 
towards  the  earth  and  the  z  axis  pointing  along  the  axis  of  rotation. 
Unfortunately,  equations  (2.2. 1-3)  do  not  hold  exactly  with  the  inertial 
reference  frame  used  by  PEP  in  this  thesis.  Since  the  lunar  moments  of 
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inertia  can  be  experimentally  determined  from  earth,  it  is  possible  to 
simultaneously  process  lunar  rotation  observations  with  satellite 
observations  such  as  those  simulated  in  this  theses.  In  this  manner, 
equations  (2.2. 1-3)  and  the  non-zero  equations  in  (2.2. 1-2)  could  be  assumed  to 
hold  exactly.  Processing  lunar  laser  corner  reflector  observations  [13,  17] 
would  provide  estimates  of  the  second  degree  coefficients  with  a  very  high 
degree  of  confidence  and  the  lunar  satellite  observations  could  be 
simultaneously  processed  for  estimates  of  the  higher  degree  coefficients. 


2.2.2  Limitations  of  the  Spherical  Harmonic  Expansion 

Although  the  spherical  harmonic  expansion  is  frequently  used  to 
describe  the  gravitational  potential  of  a  distributed  mass,  it  is  not  an  efficient 
model  for  all  uses.  Since  the  model  is  based  upon  spherical  coordinates,  it 
produces  an  uneven  resolution  of  coverage  from  the  equator  to  the  poles. 
This  uneven  resolution  can  lead  to  modeling  inefficiencies. 


Figure  2.2.2-1:  Zonal,  Tesseral,  and  Sectorial  Harmonic  Patterns 


The  spherical  harmonic  model  breaks  the  sphere  up  into  zonal, 
tesseral,  and  sectorial  patches  (Figure  2.2.2-1).  The  patches  are  separated  by 
lines  where  the  terms  are  identically  zero.  On  one  side  of  the  line,  the  term 
will  be  positive  and  on  the  other  it  will  be  negative.  The  zonal,  Jn,  terms,  are 
independent  of  longitude,  8,  and  dissect  the  globe  of  interest  into  n+1  bands 
along  lines  of  latitude.  The  tesseral  terms,  Cnm  and  &nm/  dissect  the  globe  of 
interest  into  patches  of  both  latitude  and  longitude.  A  tesseral  term  will  have 
n-m+1  sections  of  latitude  and  2m  sections  along  lines  of  longitude.  The 
sectorial  terms,  Cnn  and  Snn,  are  independent  of  latitude,  <j>,  and  dissect  the 
globe  of  interest  into  slices,  like  sections  of  an  orange. 
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Defining  a  spherical  harmonic  gravitational  model  with  a  specific 
resolution  at  the  equator  provides  an  overly  fine  longitudinal  resolution  at 
the  poles.  Obviously,  finer  surface  resolutions  require  higher  degree  models. 
The  number  of  harmonic  coefficients,  however,  grows  as  the  square  of  the 
degree,  so  increasing  the  degree  of  a  model  is  not  an  insignificant  task.  A 
more  efficient  model,  in  terms  of  coefficients  required  to  define  a  desired 
surface  resolution,  would  use  equal  size  mass  density  patches  over  the  surface 
of  the  body.  Table  2.2.2-1  lists  the  degree  and  order  of  spherical  harmonic 
models  and  the  number  of  surface  mass  density  patches  required  to  get 
various  lunar  surface  angular  resolutions  for  the  gravitational  potential.  The 
table  ignores  the  difficult  patch  layout  problem  and  resultant  inefficiencies,  so 
that  the  number  of  n°  X  n°  surface  mass  density  patches  is  the  patch  area  in 
steradians  (n2Jt2/1802)  divided  into  the  number  of  steradians  in  a  sphere  (4jc). 

Table  2. 2.2-1:  Spherical  Harmonic  and  Surface  Mass  Density  Surface  Resolution  Comparison 


Resolution 

Spherical  Harmonics 

Surface  Layer 

Degrees 

Km 

Traverse  Time 
(100  km  Alt) 

Degree  & 
Order 

#of 

Coefficients 

#  of  Patches 

22.5 

683 

442s 

16X16 

285 

82 

11.25 

341 

221s 

30X30 

957 

326 

3.297 

100 

65s 

109X109 

12,096 

3,795 

1.648 

50 

32  s 

218X218 

47,957 

15,190 

1.000 

33 

20  s 

360X360 

130,317 

41,253 

Table  2.2.2-1  lists  the  travel  time  to  cross  a  patch  of  given  size,  since  this 
factor  is  important  in  determining  the  resolution  obtainable  by  observing  a 
low-altitude  orbiter.  For  estimation  purposes,  there  should  be  at  least  two 
observations  (or  measurements)  within  the  time  it  takes  to  fly  over  a  given 
patch.  Therefore,  with  60  second  Doppler  count  intervals  of  a  satellite  in  a  100 
km  altitude  orbit,  a  spherical  harmonic  expansion  of  degree  and  order  30  is 
theoretically  possible;  more  than  three  observations  are  obtained  per  surface 
patch.  Unfortunately,  this  rule  of  thumb  does  not  account  for  smaller 
spherical  harmonic  patches  at  the  poles  or  the  lack  of  observations  during 
lunar  occultations. 

As  the  table  shows,  a  surface  mass  density  model  provides  an  economy 
in  the  number  of  coefficients  estimated,  and  might  be  a  preferred  approach  for 
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modeling  or  estimating  the  lunar  gravitational  field.  The  simulations  of  this 
thesis,  however,  were  done  with  the  spherical  harmonic  model  since  a 
surface  mass  density  model  has  not  been  validated  for  PEP.  Nonetheless, 
some  of  these  other  modeling  techniques  are  discussed  here  for  completeness. 


2.3  Localized  Surface-Layer  Gravitational  Field  Models 

Since  lunar  mascons  act  as  localized  gravitational  disturbances,  their 
modeling  requires  a  fine  degree  of  resolution  to  capture  their  high  frequency 
content.  A  very  high  degree  spherical  harmonic  model  is  required  to  model 
this  behavior,  which  requires  the  estimation  of  a  very  large  number  of 
coefficients.  As  shown  in  the  previous  section,  a  more  efficient  method  of 
modeling  this  behavior  may  be  obtained  through  models  focusing  on  the 
local,  rather  than  global,  behavior.  In  an  attempt  to  recreate  their  high 
harmonic  frequency  behavior,  mascons  have  been  modeled  by  point  masses, 
lens  shaped  mass  concentrations,  and  gravity  dipoles.  In  some  instances, 
these  mascon  models  have  been  used  to  model  the  entire  lunar  gravitational 
field,  and  in  other  cases,  they  have  been  combined  with  low  degree  spherical 
harmonic  expansions. 


2.3.1  Point  Mass  Model 

The  point  mass  model  can  be  used  to  represent  the  behavior  of  an 
individual  mascon.  The  gravitational  force  due  to  a  point  mass.  Equation 
(2.1-1),  results  from  taking  the  gradient  of  the  following  potential. 

tt/-  \  -Gm2 

U(r1)  =  jr— (23.1-1) 
lrl  r2\ 

Multiple  mascons  can  be  modeled  by  summing  the  potential  contribution  of 
each  individual  point  mass.  The  resulting  potential  model  for  n-1  point 
masses  becomes 
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Typically,  these  point  masses  are  distributed  about  known  mascon 
locations.  To  recreate  mascon  behavior,  the  location,  mass  and  depth  of  the 
point  masses  are  adjusted  to  fit  observations.  Point  masses  are  positioned 
below  the  lunar  surface  to  avoid  singular  conditions  which  would  result  as 
the  separation  distance  approaches  zero.  To  globally  model  the  moon's 
gravitational  field  with  point  masses,  the  model  should  constrain  the  total 
mass  and  center  of  mass  of  the  system  to  known  values. 


23.2  Surface  Disk  Model 


Because  point  mass  models  did  not  satisfactorily  fit  lunar  orbit 
observations,  scientists  turned  to  more  sophisticated  surface  layer 
representations.  Scientists  at  the  Aerospace  Corporation  [47]  and  the  Jet 
Propulsion  Laboratory  [4]  replaced  the  point  mass  model  with  a  surface  disk 
or  lens  shaped  model.  This  model  is  derived  from  the  potential  of  an 
ellipsoid  of  uniform  density  which  is  given  by  Equation  (2.1-4)  with  the 
boundary  of  the  ellipsoid  used  for  the  limits  of  integration.  The  density 
function  of  (2.1-4)  is  assumed  to  be  constant  and  integrates  to  die  mass  of  the 
body.  The  boundary  condition  of  an  ellipsoid  is 


(2.3.2-1) 


where  a,  b,  and  c  are  the  dimensions  of  the  ellipsoid  along  the  principal  axes 
x,  y  and  z.  The  surface  disk  or  lens  model  uses  the  specialized  case  of  an 
oblate  spheroid.  Figure  2.3.2-1  shows  a  prolate  spheroid  (a=b<c)  and  an  oblate 
spheroid  (a=b>c).  The  gravitational  attraction  for  a  point  outside  an  oblate 
spheroid  of  uniform  density  is  [34] 


y  _ 


3  Gm 


2  aV 


-e 


1+k  +  e2 
Vl +k  V  1  +k 


+  sin"1  ^ 


Vi+fc 


(2.3.2-2a) 
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where  e  is  the  spheroid's  eccentricity,  so  that  c2=a2(l-e2),  and  the  variable  k  is  a 
positive  solution  to  the  quadratic  equation 

a2k 2  +  (a2  ~  ( x 2  +  y2  +  z2  -  z2  =  0  (232-3) 


Figure  2.32-1:  Prolate  and  Oblate  Spheroids 

As  the  thickness,  c,  of  the  spheroid  approaches  zero,  the  gravitational 
attraction  of  the  disk  is  obtained.  The  resulting  forces  are  then  [47] 
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3 Gm 
a3 


(2.3.2-4) 


The  mass  of  the  disk  is  m,  a  is  its  radius,  and  the  x,  y,  and  z  coordinates  are  as 
in  Figure  2. 3. 2-1.  In  the  limit  as  the  disk's  radius,  a,  approaches  zero,  the  disk 
shaped  model  approaches  the  point  mass  model. 


The  Aerospace  Corporation  scientists  used  600  surface  disks  of  50  km 
radii  covering  the  lunar  surface  to  model  its  gravitational  field  [47].  At  JPL, 
they  used  11?  lens-shaped  mass  concentrations  placed  about  50  km  below  the 
lunar  surface  to  model  the  gravitational  field  [5].  In  the  JPL  model,  the  lens¬ 
shaped  mass  concentrations  augmented  the  moon's  central  body  attraction. 
Both  models  employed  positive  and  negative  mass  disks  in  order  to  recreate 
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the  observations  from  lunar  satellites.  Naturally,  negative  (repulsive)  mass 
surface  disks  are  not  physically  realizable,  but  they  do  model  situations  where 
a  lunar  surface  feature  is  significantly  less  dense  than  its  surroundings. 

Both  of  these  surface  disk  models  provided  better  correlation  to  lunar 
observations  than  previous  models.  Unfortunately,  they  did  not  constrain 
the  lunar  center  of  mass.  When  the  models  were  converted  to  spherical 
harmonic  models,  Ji,  Cn  and  Sn  coefficients  were  required.  Additionally, 
both  models  assumed  a  priori  knowledge  of  the  mass  concentration  locations 
and  placed  them  about  the  moon's  surface  (or  just  below  it)  on  a  grid  pattern. 
If  such  a  priori  knowledge  were  not  used,  the  estimation  process  would 
require  five  terms  to  describe  each  surface  disk:  radius  from  the  moon's 
center,  latitude,  longitude,  strength  (mass),  and  disk  radius.  A  better  model 
would  have  allowed  the  concentration's  location  to  vary  and  would  have 
constrained  the  center  of  mass.  This,  however,  would  have  involved  too 
many  variables  for  the  model  to  converge  with  the  given  lunar  orbiter  data. 


133  Gravity  Dipole  Model 

A  gravity  dipole  model  has  also  been  proposed  to  account  for  the  effect 
of  anomalous  mass  concentrations  upon  low  orbiting  bodies  [18].  A 
gravitational  dipole  consists  of  a  mass  +m  separated  from  a  fictitious  mass  -m 
by  a  distance  d  (Figure  2.3.3-1).  In  the  limit  as  d  approaches  zero,  m  is 
assumed  to  get  larger  so  that  the  product  md  remains  constant.  The  strength 
of  the  dipole,  D,  is  the  result  of  the  following  limit 

n  Um  m(r')  d  =  D(f ')  (2.3.3-1) 

Gravitational  dipoles  have  never  been  shown  to  physically  exist,  but  a 
distribution  of  these  dipoles  can  be  useful  for  modeling  the  lunar  gravity 
field. 


The  gravitational  potential  due  to  a  gravity  dipole  may  be  written  as 

D(r')  cos  6 


-ds 


(2.3.3-2) 


43 


LUNAR  GRAVITATIONAL  FIELD  ESTIMATION  AND  SATELLITE  ORBIT  PREDICTION 


where  ds  is  a  patch  of  the  surface  and  the  other  elements  are  pictured  in 
Figure  2.3.3-1. 


Gravity  dipoles  can  be  useful  modeling  tools  because  of  their  unique 
properties.  The  most  useful  property  is  that  they  produce  a  discontinuity  in 
the  tangential  component  of  the  gravitational  field  when  traversing  the 
dipole  layer.  The  gravitational  field  normal  to  the  dipole  layer,  however, 
remains  continuous  as  the  layer  is  crossed.  This  allows  gravity  dipoles  to 
model  unexplained  out  of  plane  accelerations.  Using  a  ring  of  mass  and 
gravity  dipoles  in  a  continuous  line  distribution  around  the  lunar  equator; 
Croopnick  [18]  successfully  modeled  actual  disturbing  accelerations  beneath  a 
low  altitude  orbiter.  The  dipoles  were  oriented  normal  to  the  equatorial 
plane,  so  that  they  were  used  to  account  for  out-of-plane  disturbing 
accelerations.  The  mass  ring,  meanwhile,  accounted  for  the  radial  and 
tangential  (in  plane)  components  of  the  disturbing  accelerations. 
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Chapter  Three 


Orbit  Propagation 


3.1  Numerical  Integration  Techniques 

Two  classes  of  perturbation  methods  are  used  in  celestial  mechanics  to 
determine  precise  spacecraft  orbits.  General  perturbations  generalize  the 
expressions  for  two-body  motion  to  include  disturbing  effects  of  other  bodies 
using  infinite  trigonometric  series  expansions  and  integrate  these  series  term 
by  term.  Special  peturbations  use  numerical  methods  for  deriving  the 
disturbed  orbit  by  direct  integration  of  the  rectangular  coordinates  or  a  set  of 
osculating  orbital  elements  [10]. 

Orbit  propagation  in  PEP  is  performed  using  the  numerical  methods  of 
special  perturbations.  Numerical  orbit  determination  techniques  are  prefered 
because  of  the  ease  of  implementation  and  the  accuracy  of  solutions.  The 
growing  capabilities  and  increasing  speeds  of  modem  digital  computers  have 
significantly  increased  this  method's  accuracy  and  utility.  Using  Cowell's 
method,  PEP  integrates  the  equations  of  motion  in  rectangular  coordinates 
fixed  in  inertial  space.  Section  3.2  discusses  the  units  and  coordinate  frames 
used  in  this  study's  analyses.  The  equations  of  motion  and  the  equations  for 
the  partial  derivatives  of  motion  with  respect  to  gravity  harmonic  coefficients 
and  initial  osculating  orbital  elements,  as  used  in  PEP,  are  covered  in  Sections 
3.3  and  3.4. 

Two  of  PEP'S  numerical  integration  techniques  were  used  for  this 
study:  the  Adams-Moulton  and  Nordsieck  methods.  The  Adams-Moulton 
constant  step  size  integration  technique  with  11th  differences  was  used  for  the 
propagation  of  near-circular  lunar  orbits  [15].  This  technique  uses  predictor- 
corrector  techniques  to  accurately  extrapolate  forward  in  time  a  satellite's 
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position  and  velocity.  Its  integration  step  size  should  be  smaller  then  l/100th 
of  the  satellite's  orbital  period  to  ensure  numerical  stability.  An  even  smaller 
step  size  must  be  used  to  accurately  sample  higher  degree  gravitational 
harmonics.  A  very  conservative  rule  of  thumb  is  to  use  a  step  size  of  l/(100n) 
of  an  orbital  period  to  simulate  an  n  X  n  degree  and  order  spherical  harmonic 
model.  In  PEP  the  step  size  is  2m  days,  where  m  is  a  negative  number.  A  two 
hour  lunar  orbit  using  an  m  =  -10  step  size  will  have  over  170  steps  per  orbit, 
ensuring  numerical  stability.  An  m  =  -14  step  size  for  the  same  orbit  will 
have  over  1,365  steps  per  orbit,  which  should  adequately  model  a  13  X  13 
spherical  harmonic  gravity  model. 

Elliptical  orbits  were  propagated  by  the  Nordsieck  fifth  degree  variable 
step  size  technique  [38].  Since  the  technique  is  self  starting,  it  is  used  to  start 
the  Adams-Moulton  method.  This  technique  predicts  the  orbit  ahead  using  a 
fifth  degree  polynomial  whose  coefficients  are  approximations  to  derivatives 
of  the  function  being  integrated.  The  integration  output  file  therefore 
contains  the  satellite  position,  velocity,  acceleration,  and  jerk  [8].  The  variable 
step  size  uses  a  smaller  step  near  periapse,  where  quantities  change  more 
rapidly.  For  highly  elliptical  orbits,  this  is  a  more  efficient  integration 
technique,  since  constant  step  size  methods  propagate  the  entire  orbit  with 
the  smallest  required  step  size. 

Despite  the  integration  step  size  used,  integration  quantities  are  written 
to  an  output  file  using  a  different  step  size.  This  output  step  size  is  generally 
two  times  the  integration  step  size.  For  low  frequency  orbital  disturbances, 
less  frequent  output  step  sizes  can  be  used.  A  satellite's  position  and/or 
velocity  are  determined  at  specific  observation  times  by  interpolating  from 
the  satellite's  integration  output  file.  Everett  eighth  difference  interpolation 
is  performed  on  constant  step  size  integration  files.  This  method  fits  a  ninth 
degree  polynomial  through  the  ten  output  times  surrounding  the 
observation  time.  Hermite  interpolation  is  performed  on  variable  step  size 
integration  files.  This  method  uses  a  fifth  degree  polynomial  agreeing  with 
position,  velocity,  and  acceleration  at  the  two  output  times  surrounding  the 
observation  time  [8]. 
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3.2  Units  and  Coordinate  Systems 

The  numerical  integration  methods  of  special  perturbations  require 
that  units  of  time  and  the  gravitational  constant  are  precisely  defined.  The 
time  unit  in  PEP  is  the  coordinate  time  day  defined  in  terms  of  atomic  time  in 
Appendix  B.  The  distance  unit  is  the  Astronomical  Unit  (AU),  defined  by 
setting  tiie  square  root  of  the  gravitational  constant  times  the  mass  of  the  sun 
to  the  Gaussian  value  (see  Appendix  B). 

Numerical  integration  of  the  satellite  equations  of  motion  are 
performed  in  an  inertial  Cartesian  coordinate  system.  PEP'S  current  inertial 
coordinate  system  is  based  upon  the  mean  equinox  and  equator  of  the  earth  of 
1950.0,  or  Julian  Date  2,433,282.423.  This  coordinate  system  uses  the  earth's 
rotation  axis  at  this  time  as  the  z  axis,  with  the  x  axis  along  the  1950.0  mean 
equinox  pointing  towards  the  constellation  Aries,  and  the  y  axis  completing  a 
right-handed  coordinate  system.  The  transformations  between  this 
coordinate  system  and  those  fixed  in  the  moon  and  earth  are  discussed  in 
Appendix  C. 

For  lunar  satellite  propagations,  the  origin  of  this  coordinate  system  is 
placed  at  the  moon's  center  of  mass.  The  coordinates  of  perturbing  bodies 
during  an  integration  are  determined  from  the  coordinates  of  the  earth-moon 
barycenter  relative  to  the  sun,  of  the  moon  relative  to  the  earth,  and  of 
planets  relative  to  the  sun  calculated  by  interpolation  from  an  n-body  file 
supplied  from  the  Smithsonian  Astrophysical  Observatory  (SAO).  This  n- 
body  file  is  based  upon  the  SAG'S  fit  to  observational  data.  The  moon's 
coordinates  are  based  upon  formulas  from  Brown's  lunar  theory. 

To  analyze  the  PEP  propagated  orbits,  auxiliary  software  was  written  to 
transform  the  inertial  integration  Cartesian  coordinates  into  a  selenographic 
coordinate  system.  The  selenographic  coordinate  system  is  also  centered  at 
the  lunar  center  of  mass,  but  its  z  axis  is  the  lunar  rotation  axis,  the  x  axis 
points  toward  the  earth,  and  the  y  axis  completes  the  right  hand  coordinate 
system  (see  Appendix  C.l).  One  auxiliary  software  program  computes 
selenographic  orbital  elements  and  satellite  ground  tracks  from  an  inertial 
integration  file.  Initial  conditions  for  lunar  orbits  were  chosen  in  the 
selenographic  coordinate  system.  The  orbital  angles  (i,  Q,  co)  were  then 
transformed  to  the  inertial  integration  coordinate  frame  in  a  second  auxiliary 
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software  program  using  the  transpose  of  the  transformation  matrix  used  in 
the  first  program. 


For  this  thesis,  inertial  Cartesian  coordinate  frame  positions  will  be 
referred  to  as  X,  and  in  selenographic  Cartesian  coordinates  as  U,  where  the 
two  are  related  by  the  transformation  matrix  R(t)  as  follows  (Appendix  C.l) 
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Additionally,  the  selenographic  Cartesian  coordinates  are  related  to  the 
selenographic  spherical  coordinates  (r,$,8)  by  the  relationship 


U 

r  cos  6  cos  (p 

V 

— 

r  sin  0  cos 
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so  that  the  inertial  coordinate  frame  can  be  related  to  the  selenographic 
spherical  coordinates  used  in  the  spherical  harmonic  expansion  by  the 
relations 


r  =  |x|  =  (XT£)*  =  |Q|  02-4) 

sin0  =  cos0  =  ^l-sin2  <p  (32-5) 

°>s0  =  7s!m(RiTX)  sine  =  7^y(RjX)  (3.2-6) 

These  transformations  are  used  in  PEP's  internal  transformation 
routines,  in  the  two  auxiliary  software  programs,  and  in  the  partial 
derivatives  of  satellite  motion  equations  covered  in  Section  3.4. 
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3.3  Equations  of  Satellite  Motion 


Because  of  the  nature  of  the  equations  of  motion  for  a  satellite,  PEP 
integrates  an  augmented  state  vector, 

x;,  =  f  §W-  ]  (3.3-1) 


x;=KJ 

where  Xw  is  the  position  of  body  b  (satellite)  with  respect  to  body  1  (the  moon) 
defined  by  (3.2-1)  and  Xw  is  its  time  derivative: 

<>  _  dXbl  _  r .  -|T  /o 


4=^  =  [*  y  if 


(3.3-2) 


For  the  following  equations,  X,  X,  and  X*  will  all  be  used,  although 
the  reader  should  realize  that  "state"  refers  to  the  augmented  state  vector,  X*. 

The  equations  of  motion  for  a  lunar  satellite  relative  to  the  moon  may 
be  written  as: 


<*2XW  _  Xw 


=  — GA lj  — j — h  ffj  +  'F  +  He  +  other  forces 
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(3.3-3) 


with  the  initial  condition  at  t=to 


Xw  -  xwo 


(33-4) 


In  this  thesis,  the  "other  forces"  (radiation  pressure,  gas  leaking, 
thruster  firing)  are  ignored,  but  can  be  included  when  fitting  to  real 
observations.  In  Equation  (3.3-3),  the  H /  and  H*  terms  are  the  effects  of 
gravitational  harmonics  (zonal  and  tesseral)  for  the  moon  and  earth 
respectively,  and  the  'V  term  is  the  point  mass  perturbing  accelerations  of 
other  bodies  upon  the  spacecraft  and  the  subscripts  e,  s,  and  p  refer  to  the 
earth,  sun  and  planets. 


(3.3-5) 


For  the  simulations  run,  the  earth  and  sun  perturbing  attractions  were 
included.  The  effect  of  all  other  bodies  was  considered  negligible  for  these 
simulations. 


3.3.1  Effect  of  Lunar  Gravitational  Harmonics 


The  acceleration  on  a  satellite  due  to  the  gravitational  potential  U  of 
the  moon  is  obtained  by  taking  the  gradient  of  the  potential  (2.1-3).  Since 
Equation  (3.3-3)  already  accounts  for  the  central  body  attraction,  the  effect  of 
the  remaining  harmonics  is  derived  from  (2.2-9)  without  the  central  body 
term.  The  resulting  acceleration  is 
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where  X  is  the  ratio  of  pz/r m,  Nz  is  the  degree  of  zonals  used  in  the  spherical 
harmonic  expansion  and  N*  is  the  degree  and  order  of  tesseral  harmonics  in 
the  expansion.  The  recursive  formulas  for  Pn,  P'n,  Pnm,  and  P'm  in  terms  of 
the  argument  sinij)  which  have  been  coded  in  PEP  are  given  in  Appendix  A. 

The  partial  derivatives  of  the  selenographic  spherical  angles  in  (3.3. 1-1) 
are  obtained  by  differentiating  Equations  (3.2-5)  and  (3.2-6).  This 
differentiation  results  in  the  following  relationships 
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Equations  (3.3.1-3)  and  (3.3.1-4)  can  be  simplified  by  multiplying  the 
first  by  -sin0,  the  second  by  cos0,  and  then  adding  the  two  equations.  This 
yields  the  following  formula  for  the  partial  derivative  with  respect  to  0. 
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(3.3.1-5) 


Equation  (3.3.1-1)  is  singular  at  (j>  =  ±90°  for  PEP'S  algorithms.  This  is 
generally  not  a  problem  since  it  is  highly  unlikely  that  a  satellite  numerical 
integration  step  would  land  exactly  over  a  pole.  One  simulation,  however, 
used  a  polar  location  as  an  initial  condition  and  PEP  was  unable  to  propagate 
its  orbit.  The  initial  mean  anomaly  for  this  orbit  was  altered  by  one  degree 
and  the  integration  proceeded  without  any  further  difficulties. 


3.3.2  Effect  of  Other  Gravitational  Body  Harmonics 

PEP  has  the  capability  to  include  other  (non-central)  gravitational  body 
perturbing  effects  upon  a  satellite's  motion.  This  enables  PEP  to  handle 
spacecraft  fly-by  missions.  PEP  can  include  the  effects  of  a  destination,  or 
target,  body's  harmonics  upon  the  motion  of  a  satellite  traveling  towards  one 
body  but  which  is  within  the  sphere  of  influence  of  another  body.  Using  this 
feature,  PEP  can  also  include  the  earth's  J2  harmonic  upon  a  lunar  orbiting 
satellite.  Higher  order  terms  can  also  be  included  for  integration  accuracy,  but 
are  rarely  required. 

The  effect  of  earth  perturbations  on  the  motion  of  a  lunar  satellite 
relative  to  the  moon  (H*  from  (3.3-3))  is  the  difference  between  the  effect  of 
the  earth  on  the  satellite  (Hj,e)  and  the  effect  of  the  earth  on  the  moon  (H ie). 
Hfc*  is  given  by  Equation  (3.3.1-1)  with  a  change  of  subscripts,  replacing 
subscript  1  with  e.  The  effect  of  earth  harmonics  upon  the  moor  is  calculated 
by  replacing  the  GMj  terms  in  (3.3.1-1)  by  GMc,  (the  subscript  c  refers  to  the 
earth-moon  barycenter),  the  subscript  1  by  e,  and  the  subscript  b  by  1.  The  H* 
term  in  (3.3-3)  is  then  given  by 

He  =  Hte-Hfc  (3.3.2-1) 
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Since  low  altitude  lunar  orbits  and  eccentric  orbits  with  apolune  on  the 
far  side  of  the  moon  were  simulated,  no  earth  perturbing  effects  were 
included.  When  processing  real  observations  to  estimate  the  lunar 
gravitational  field,  earth's  J2  harmonic  should  be  considered  since  it  is  three 
orders  of  magnitude  larger  than  the  other  earth  harmonics. 

3.3.3  Effect  of  Mass  Concentrations 

For  this  thesis,  PEP  was  modified  to  include  the  effects  of  mascons 
upon  lunar  orbiting  satellites.  Since  these  modifications  have  not  been 
integrated  into  the  SAOXs  version  of  the  software,  the  Draper  Laboratory 
modified  version  will  be  referred  to  as  PEP-D.  PEP-D  implements  mascons  by 
modeling  them  as  point  masses  (Section  2.3.1). 

Due  to  mascon  model  implementation,  the  equations  of  motion  for  a 
lunar  satellite  relative  to  the  moon  (3.3-3)  become,  in  PEP-D, 

+  +"¥  +  Ht+K  (3.3.3-1) 

at  rbl 

with  the  same  initial  conditions,  (3.3-4).  In  Equation  (3.3.3-1),  the  K  term  is 
the  effect  of  all  of  the  mascons  and  u  is  their  total  mass  fraction. 


The  mass  of  each  mascon  and  its  selenographic  position  in  spherical 
coordinates  is  entered  into  PEP-D  as  the  program  is  initialized.  The  mass  or 
strength  of  a  mascon  is  input  as  a  fraction  of  the  total  central  body  mass  and 
both  positive  and  negative  values  are  allowed.  If  Njt  is  the  total  number  of 
mascons,  the  total  mass  fraction  of  the  mascons  is  then  calculated  by 


V 


0<  U<1 


(3.3.3-2) 


The  original  central  body  term  from  (3.3-3)  is  reduced  by  the  factor  (1-u) 
to  conserve  mass  in  the  lunar  system.  This  will  allow  the  lunar  mascon 
model  to  behave  identically  to  the  central  force  model  far  from  the  moon. 
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The  acceleration  on  a  satellite  due  to  the  gravitational  potential  of  N* 
mascons  is  obtained  by  taking  the  gradient  of  the  point  mass  potential  (2.3.1- 
2).  The  resulting  acceleration  is 

Nk  Y 

K  =  0.3.3-3) 

i=l  Tbi 

From  the  input  spherical  mascon  selenographic  coordinates,  PEP-D 
determines  their  inertial  Cartesian  coordinates  Xa  using  Equation  (3.2-3)  and 
the  inverse  of  Equation  (3.2-2).  The  coordinates  of  the  lunar  satellite  relative 
to  mascon  i  and  their  separation  are  then 

Xbi  —  Xbi  —  Xu  (3.3.3-4) 

rw=|Xw|  (3.3.3-5) 

Unfortunately  the  mascon  implementation  in  PEP-D  does  not 
constrain  the  lunar  center  of  mass.  As  mentioned  in  Section  2.3.2,  when 
converting  a  surface  layer  model  with  an  unconstrained  center  of  mass  to  a 
spherical  harmonic  expansion,  first  degree  harmonic  coefficients  Ji,  Cn,  and 
Sn  need  to  be  determined.  Rather  than  modify  PEP-D  to  estimate  first  degree 
coefficients,  this  thesis  uses  mascon  models  in  which  the  lunar  center  of  mass 
is  not  disturbed.  If  any  mascons  are  used,  PEP-D  calculates  the  central  body 
center  of  mass  by  the  formula 

Nk 

Xcm  =XmA  (3.33-6) 

i= 1 

When  the  somewhat  arbitrary  mascon  "truth"  model  was  created  (Section 
6.4),  this  feature  was  used  to  determine  the  placement  and  mass  of 
"balancing"  mascons  to  preserve  the  lunar  center  of  mass. 


3.4  Partial  Derivatives  of  the  Satellite  Motion  Differential  Equation 

The  partial  derivative  of  Equation  (3.3-3)  with  respect  to  a  parameter  a 
yields  the  variational  equations  that  are  numerically  integrated  in  PEP  along 
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with  the  equations  of  motion.  The  parameter  a  could  be  one  of  the  orbit's 
initial  conditions,  the  harmonic  coefficients  of  the  gravitational  field,  or  other 
parameters  of  interest.  Taking  into  account  the  effects  of  mascons,  the  partial 
derivative  of  Equation  (3.3.3-1)  with  respect  to  the  same  parameter  a  yields 
PEP-D's  variational  equations 
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(3.4-1) 


with  the  initial  condition  at  t  =  to 
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(3.4-2) 


Unless  estimating  the  mass  of  the  moon,  the  first  term  in  (3.4-1)  is  zero. 
Additionally,  since  the  mascon  masses  were  not  estimated  in  this  thesis,  the 
ratio  u  was  considered  constant. 


The  effect  of  perturbing  body  attractions  (point  mass  approximation) 
upon  the  partial  derivative  with  respect  to  a  is  obtained  by  differentiating 
Equation  (3.3-5),  yielding 
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(3.4-3) 


3.4.1  Effect  of  Lunar  Spherical  Harmonics  on  the  Partial  Derivatives 

The  effect  of  the  lunar  gravitational  harmonics  upon  the  partial 
derivatives  of  satellite  motion  with  respect  to  a  parameter  a  is  obtained  by 
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differentiating  Equation  (3.3.1-1)  term  by  term.  After  differentiating,  like 
series  summations  can  be  regrouped  for  numerical  algorithms.  The  resulting 
equation  takes  the  form 


dHl_  ft,  d(GM,) 
da'GM,  da  rj  ^ 


(3.4.1-1) 


where  2  is  the  series  expansion  for  the  zonal  terms  and  ?  is  the  series 
expansion  for  the  tesseral  terms.  The  first  term  in  (3.4.1-1)  is  set  to  zero  unless 
attempting  to  estimate  the  mass  of  the  moon.  2,  the  zonal  series  expansion, 
is  provided  by  the  following  expression: 
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(3.4.1-2) 

The  first  term  in  the  above  equation  is  zero  since  the  zonal  coefficients  are 
constant. 

After  grouping  like  terms,  the  series  expansion  for  the  tesserals,  f, 
reduces  to  the  following  expression: 
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(3.4.1-3) 


The  first  two  terms  in  this  expansion  are  zero  since  the  tesseral  spherical 
harmonic  coefficients  are  constant.  The  expansion  for  Ti  and  T2  are  then 
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The  recursive  formulas  for  calculating  Pn/  P'„,  P",  Pnm,  P'm,  and  P"m 
in  terms  of  their  argument  sin$  are  included  in  Appendix  A. 

Additionally,  based  upon  the  transformation  between  selenographic 
spherical  coordinates  and  the  inertial  frame  (3.2-4),  (3.2-5)  and  (3.2-6),  the 
following  partial  derivatives  are  obtained. 
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PEP  has  the  option  of  using  different  degrees  of  spherical  harmonic 
models  for  calculating  the  partial  derivatives  in  Equation  (3.4-1)  than  for 
calculating  the  satellite  motion  in  Equation  (3.3-3)  or  (3.3.3-1).  Since 
calculating  the  effect  of  the  harmonics  is  such  a  tedious  process,  calculating 
the  partial  derivative's  harmonic  effect  with  a  lower  degree  spherical 
harmonic  model  can  save  quite  a  bit  of  computer  time.  This  option  works 
well  for  high  altitude  earth  satellites  where  J2  dominates  the  other  harmonic 
effects.  This  option  was  found  impractical  for  low  altitude  lunar  satellites,  so 
the  same  degrees  Nz  and  Nt  were  used  as  the  summation  upper  limits  in 
Equations  (3.4-1)  as  in  (3.3.3-1). 


3.4J2  Effect  of  Earth  Spherical  Harmonics  on  the  Partial  Derivatives 

The  effect  of  earth  harmonics  on  the  partial  derivatives  of  satellite 
motion  with  respect  to  a  parameter  a  is  calculated  by  differentiating  Equation 
(3.3.2-1).  This  calculation  involves  differentiating  two  factors  of  the  form  of 
Equation  (3.4.1-1)  since  the  effect  of  the  earth  harmonics  on  the  moon  must  be 
subtracted  from  the  effect  of  the  earth  harmonics  upon  the  lunar  satellite.  For 
this  thesis,  the  effect  of  earth  harmonics  upon  the  partial  derivatives  were 
neglected  in  Equation  (3.4-1)  because  all  of  the  orbits  propagated  were  low 
altitude  ones: 

9He  -  / 
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3.13  Effect  of  Mass  Concentrations  upon  the  Partial  Derivatives 


The  effect  of  mass  concentrations  on  the  partial  derivatives  of  satellite 
motion  with  respect  to  a  parameter  a  is  obtained  by  differentiating  Equation 
(33.3-3)  assuming  that  the  total  lunar  mass  remains  constant. 
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The  above  equation  was  not  coded  into  PEP-D.  The  partial  equations 
were  calculated  on  orbit  fitting  runs  and  all  of  the  orbit  fitting  was  done  to 
spherical  harmonic  models.  The  spherical  harmonic  plus  mascon  truth 
model  was  run  to  generate  the  satellite  observations  for  each  gravitational 
sensing  method,  for  which  partial  derivatives  were  not  required.  If  PEP-D 
was  used  to  estimate  coefficients  in  a  spherical  harmonic  plus  mascon  model, 
this  equation  would  be  required  and  could  be  used  in  an  estimation  fit  which 
varied  the  mascon  strength  and  location. 


3.4.4  Effects  of  Initial  Conditions  on  the  Partial  Derivatives 

The  partial  derivatives  of  the  satellite's  initial  conditions  with  respect 
to  a  from  Equation  (3.3-4)  (used  as  initial  conditions  (3.4-2))  are  zero  unless 
the  parameter  itself  is  a  satellite  initial  condition.  For  the  cases  where  a  is  a 
Cartesian  coordinate  initial  condition,  the  partial  derivative  is  obtained  from 
the  following  relation. 


dXbip  _  T  _ 

ax'm 

Rather  than  solve  the  partial  derivatives  with  respect  to  inertial 
Cartesian  initial  conditions,  PEP  uses  the  partial  derivatives  of  the  initial 
osculating  elliptical  orbital  elements  (ao,  eo,  io,  Do,  coo.  Mo)  or  (ao,  eo,  io,  Do, 
(D+g))o,  (D+cd+M)o).  The  osculating  orbit  is  the  one  which  would  result  if  the 
disturbing  forces  in  Equations  (3.3-3)  or  (33.3-1)  were  instantly  "turned  off" 
and  the  satellite's  motion  continued  along  the  two-body  orbit  defined  by  the 


(3.4.4-1) 
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central  body  force.  These  orbital  elements,  like  the  Cartesian  coordinate  state 
vector,  completely  describe  the  satellite's  initial  state.  Calculating  the  partial 
derivatives  with  respect  to  the  initial  osculating  elliptical  orbital  elements  is 
numerically  more  efficient  because  only  the  partial  derivative  with  respect  to 
the  initial  osculating  semi-major  axis  grows  with  time.  Reference  [6]  contains 
the  formulas  for  the  partial  derivatives  of  the  initial  conditions  with  respect 
to  the  initial  osculating  elliptic  orbital  elements  in  the  integration  coordinate 
frame. 


3.4.5  Checking  Partial  Derivatives  by  the  Difference  Method 


To  verify  that  partial  derivatives  were  being  calculated  correctly,  they 
were  checked  using  a  finite  difference  method.  This  method  was  used  to 
check  both  position,  X,  and  observation  partials.  For  these  checks,  orbits 
and/or  observations  with  partials  were  generated  with  two  different  values  of 
the  parameter  a  (ai  and  (X2).  The  finite  difference  method  then  verifies  the 
partial  derivatives  of  the  quantity  of  interest,  T  (scalar  position  coordinate  or 
theoretical  observation),  by  verifying  the  following  equality. 
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This  check  allows  coding  errors  to  be  detected  and  verifies  the  proper 
execution  of  the  program.  Once  coding  errors  are  fixed  and  the  program  is 
being  run  properly  (appropriate  control  inputs  are  used),  the  orbit  fitting 
process  generally  converges  to  a  solution.  If  the  program  is  running  properly 
but  has  difficulty  converging  upon  a  solution,  this  indicates  that  parameters 
are  too  highly  correlated,  the  observations  do  not  provide  adequate 
observability,  or  the  fit  model  is  an  inadequate  representation  of  the  observed 
behavior. 
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Chapter  Four 


Satellite  Observations 


4.1  Observation  Simulation  and  Processing 

Knowledge  of  the  world  in  which  we  live  is  based  upon  observations, 
or  measurements,  of  physical  phenomenon.  Theories  that  attempt  to  explain 
and  model  natural  occurrences  are  developed  to  explain  previous 
observations  and  are  validated  by  accurately  predicting  future  behavior.  The 
discovery  of  die  planet  Neptune  is  an  example  of  the  impact  observations 
have  upon  the  development  and  refinement  of  theoretical  models.  After  the 
discovery  of  Uranus,  astronomers  could  not  reconcile  observations  of  the 
planet's  motion  with  theoretical  predictions.  English  astronomer  John  Couch 
Adams  and  French  astronomer  Urbain-Jean-Joseph  Le  Verrier,  independently 
studying  the  motion,  concluded  that  Uranus'  behavior  was  due  to  a  planet 
beyond  it.  Using  a  new  model,  both  astronomers  were  able  to  predict  die 
location  of  this  ultra-Uranian  planet,  and  shortiy  thereafter  Neptune  was 
discovered  [10,  pp.  472-473]. 

The  Planetary  Ephemeris  Program  (PEP)  was  designed  to  process 
observations  of  the  sun,  moon,  planets,  stars,  and  spacecraft  and  to  further 
scientific  knowledge  through  this  observation  processing  [8].  PEP  processes 
observations  of  interplanetary,  earth,  and  lunar  spacecraft  using  observation 
files.  These  files,  based  on  actual  astronomical  observations  or  simulated 
observations,  are  made  up  of  observation  series.  Observation  series  are 
limited  to  a  single  sensing  method  (Section  4.2)  and  a  pair  of  observation 
types  (Section  4.3).  Multiple  observation  methods,  types,  periods,  and 
frequencies  are  accommodated  within  PEP  by  either  including  more  than  one 
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observation  series  in  an  observation  file  or  by  using  more  than  one 
observation  file  at  a  time. 

To  optimally  estimate  a  spherical  harmonic  expansion  fit  to  the  lunar 
gravity  field  for  this  thesis,  PEP  generated  simulated  "truth"  model 
observations  based  upon  the  observation  geometry  and  measurement  type 
being  evaluated.  Observation  series  for  each  observation  geometry  and 
measurement  type  were  created  using  PEP's  dummy  observation  feature 
where  output  observation  files  were  created  from  knowledge  of  earth-based 
observing  sites,  and  observing  and  observed  satellite  states  determined  from 
numerical  orbit  propagations  which  used  the  gravitational  "truth"  model. 
These  observation  series  included  each  method's  associated  measurement 
error  -  an  indication  of  the  statistical  accuracy  of  any  given  measurement. 
The  error  values  used  were  based  upon  published  or  calculated  instrument 
capabilities.  During  the  fitting  process,  the  observation  files  were  read  and 
these  errors  were  used  to  weight  the  observations. 

Treating  the  "truth"  model  observations  as  real,  PEP  fit  a  spherical 
harmonic  model  to  these  observations.  Using  first  guesses  for  the  fit  model 
coefficients  and  satellite  initial  osculating  orbital  elements,  the  equations  of 
motion  and  the  equations  for  the  partial  derivatives  of  the  motion  with 
respect  to  these  quantities  were  numerically  integrated.  From  the  numerical 
integrations,  the  theoretical  observation  values  were  determined. 
Observation  residuals,  the  difference  between  the  "true"  observations  and  the 
theoretical  values  calculated,  and  the  observation  partial  derivatives  were 
computed  and  written  to  an  observation  output  file.  This  file  was  later  used 
to  update  the  guess  for  the  gravitational  coefficients  and  orbital  initial 
conditions. 

There  was  no  measurement  noise  added  to  the  "truth"  model 
observations  in  this  thesis.  White  measurement  noise  could  be  added  using  a 
random  number  generator,  »/ith  any  non-zero  measurement  biases  estimated 
by  PEP. 

Both  "truth"  model  and  theoretical  observation  values  depend  upon 
the  location  of  the  observing  site,  the  observation  time  (signal  receive  time), 
the  location  of  the  observed  body  at  reflection  time,  and  any  distortion  factors. 
PEP  determines  the  receiving  site's  position  or  state  at  signal  receive  time  in 
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the  inertial  Cartesian  coordinate  frame  based  upon  integration  or  n-body  files. 
For  observation  processing  in  ds-lunar  space,  position  is  determined  relative 
to  the  center  of  the  earth.  For  extremely  precise  observations,  PEP  can 
perform  the  light  time  iterations  with  the  center  of  mass  of  the  solar  system 
as  the  center  of  the  inertial  coordinate  system. 

The  distance  between  the  observed  body  at  a  guessed  reflection, 
transpond,  or  transmission  time  and  the  observing  site  at  receive  time  is 
converted  to  light  time,  and  the  reflection,  transpond,  or  transmission  time  is 
adjusted  until  the  difference  between  the  reflection,  transpond,  or 
transmission  time  and  receive  time  is  sufficiently  close  to  the  distance 
between  the  bodies  converted  to  light  seconds. 

For  two  way  observation  signals,  this  process  is  performed  twice.  The 
first  iteration  is  used  to  determine  the  reflection  time  and  state  of  the 
observed  body.  The  second  step  uses  the  reflection  time  (less  any  systemic 
delays)  and  state  of  the  observed  body  to  compute  the  signal  source's  send 
time  and  position.  Bent  pipe  observations  perform  these  light  time  iterations 
as  many  times  as  necessary  to  recreate  the  path  of  the  electromagnetic  signals 
from  signal  receive  time  back  to  its  origin. 


4.2  Observation  Methods 

Astronomical  observations  are  collected  using  several  different 
methods.  The  most  common  observations  are  taken  from  earth  sites.  These 
can  be  observations  of  celestial  bodies  or  man-made  spacecraft  in  earth,  lunar, 
or  interplanetary  orbits.  Additionally,  satellite  based  observations  are  made  of 
other  satellites  or  of  sites  on  the  earth  or  lunar  surface.  Finally,  bent  pipe 
observations  are  made  in  which  a  signal  is  passed  among  multiple  satellites 
and  sites.  PEP  can  simulate  observations  using  each  of  these  methods  and  can 
also  model  physical  effects  which  affect  these  observations. 
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42.1  Earth-Based  Observations 

PEP  contains  the  coordinates  of  several  earth  based  sites  to  process 
astronomical  observations.  DSN  sites  are  included  in  cylindrical  coordinates 
and  most  other  locations  are  in  spherical  coordinates.  Observation  sites  not 
already  included  can  be  added  to  PEP.  This  will  permit  the  creation  of  lunar 
observation  sites.  If  a  lunar  base  or  navigation  site  were  established,  PEP 
could  process  lunar-based  observations  in  the  same  manner  used  to  process 
earth-based  observations,  after  some  straightforward  software  modifications. 

When  an  earth-based  site  is  involved  in  an  observation,  the  site's 
inertial  location  with  respect  to  the  center  of  the  earth  must  be  determined  to 
calculate  the  theoretical  value  of  the  observations.  Given  an  observation 
receive  or  send  time  in  UTC  time  (as  disseminated  by  the  U.S.  Naval 
Observatory  Time  Service  Radio  Station,  WWV,  and  other  time  services), 
PEP  determines  the  coordinate  time  (CT),  UT1  time,  and  earth  wobble 
coordinates  from  look-up  tables.  International  Astronomical  Union  (IAU) 
expressions  for  the  sidereal  time  as  a  function  of  UT1  time  and  the  earth 
precession-nutation  matrices  as  functions  of  CT  time  are  then  evaluated.  PEP 
transforms  the  observing  site's  earth  fixed  coordinates  to  the  integration 
frame  coordinate  system  using  the  transformations  described  in  Appendix 
C.2.  PEP  can  also  calculate  the  partial  derivatives  of  the  observing  site's 
integration  frame  coordinates  with  respect  to  cylindrical  or  spherical 
coordinates.  This  feature  allows  PEP  to  improve  the  estimate  of  the 
observation  site's  location  in  earth  fixed  coordinates  as  part  of  the  process.  As 
lunar  sites  are  established,  PEP  can  therefore  survey  the  site's  location  by 
processing  satellite-,  earth-,  and  lunar-based  observations. 


4.2.2  Satellite-Based  Observations 

PEP  also  allows  observations  when  a  satellite  is  the  signal  receiving 
and/or  sending  site.  In  this  case,  PEP  uses  the  satellite's  integration  file  to 
determine  its  state  and  the  partial  derivatives  of  the  state  at  receive, 
bounce /transpond,  or  send  time.  When  using  a  lunar  satellite  for 
observations,  the  satellite's  state  and  partial  derivatives  are  translated  from 
lunar-  to  earth-centered  coordinates.  The  lunar  satellite's  state  and  any  partial 
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derivatives  at  receive,  bounce/ transpond,  or  send  time  are  determined  using 
the  interpolation  methods  described  in  Section  3.1. 


4L2.3  Bent  Pipe  Observations 

PEP  also  processes  observations  in  which  signals  pass  among  several 
satellites  and  sites.  This  coding  has  been  used  to  recreate  two  different 
ranging  observations  [8].  In  the  first  method,  an  earth  site  transmitted  a 
ranging  signal  to  an  earth  satellite  which  then  transponded  the  signal  to  a 
second  earth  satellite.  This  second  satellite  then  transponded  the  signal  to  an 
earth  receiving  site.  The  second  bent  pipe  method  used  the  same  signal  path 
from  the  originating  earth  site  through  the  second  satellite,  but  then  the 
signal  was  transponded  back  to  the  first  satellite  and  then  transponded  to  an 
earth  receiving  site.  In  both  cases  the  earth  sending  and  receiving  sites  can  be 
the  same  or  different  sites.  PEP-D  can  process  these  same  bent  pipe 
observations  methods  for  lunar  orbiting  spacecraft  with  earth-based  sending 
and  receiving  sites.  Unfortunately  these  bent  pipe  methods  currently  only 
process  ranging  observations  (Section  4.3.6).  To  accommodate  bent  pipe  range 
rate  measurements,  PEP  will  have  to  be  modified. 

A  bent  pipe  method  using  two  two-way  coherent  links  has  been 
proposed  to  observe  satellite  motion  on  the  far  side  of  the  moon  (Section  1.4). 
For  this  method,  an  earth-based  observation  site  sends  a  high  frequency  signal 
to  a  lunar  satellite.  This  Doppler  shifted  received  frequency  is  used  to 
generate  a  medium  frequency  signal  for  a  two-way  coherent  Doppler  loop 
between  two  lunar  satellites.  The  returned  Doppler  shifted  medium 
frequency  is  then  modulated  onto  the  received  high  frequency  signal  and 
transmitted  to  earth,  completing  a  second  coherent  loop  [40].  Since  this 
observation  method  was  not  available  for  range  rate  measurements,  it  was 
simulated  by  using  two  separate,  unrelated  coherent  Doppler  loops.  Since  the 
proposed  method  depends  on  a  single  frequency  source,  the  simulated 
method  used  a  perfect,  non-drifting  frequency  source  for  the  satellite-to- 
satellite  loop.  This  method  obtained  range  rate  observables  equivalent  to  the 
proposed  method. 
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4.2.4  Observation  Interruptions 

All  observations  involve  the  transmission  of  an  electromagnetic  signal 
(light  or  radio  waves)  between  sending  and  receiving  sites  and  satellites. 
Observations  therefore  depend  upon  a  clear  line-of-sight  between  the  sites 
and  satellites.  Multiple  signal  path  methods  depend  upon  a  clear  line-of-sight 
for  each  path.  For  earth  observations  of  a  lunar  satellite,  the  line-of-sight  can 
be  interrupted  by  either  the  earth  or  the  moon.  The  first  case  occurs  if  the 
line-of-sight  passes  below  the  horizon.  In  the  second  case,  the  moon  occults 
the  line-of-sight  when  the  lunar  satellite  passes  behind  the  moon.  PEP 
models  both  of  these  observation  interruptions  and  deletes  observations  from 
dummy  "truth"  model  observation  series. 

Given  the  radius  of  the  occulting  body,  PEP  determines  whether  the 
satellite  is  occulted  by  its  central  body  using  vector  descriptions  of  the 
observing  site,  observed  satellite,  and  central  body  locations.  When  the 
observed  body  becomes  occulted  by  its  central  body,  observations  cease  until 
the  observed  body  returns  to  view.  For  the  gravitational  sensing  methods 
simulated,  this  feature  was  used  with  a  lunar  radius  slightly  larger  than  the 
lunar  radius  to  eliminate  poor  quality  observations  and  account  for  satellite 
acquisition  difficulties  at  the  edges. 

Given  a  limiting  elevation  angle,  as  depicted  in  Figure  4.3.1-2,  PEP 
determines  whether  the  line-of-sight  has  passed  below  the  observing  site's 
horizon  using  vector  descriptions  of  the  line-of-sight  and  the  vector  normal 
to  the  observing  site.  PEP  discontinues  observations  while  the  line-of-sight  is 
below  the  required  elevation  angle.  This  feature  was  not  used  in  the  lunar 
gravity  model  estimation  simulations  since  there  is  at  least  one  observing  site 
facing  the  moon  at  all  times  (these  sites  are  spaced  around  the  earth's 
longitude).  For  this  thesis,  the  simulations  used  a  single  observing  site  which 
made  observations  through  a  "transparent"  earth  rather  than  vising  several 
observing  sites  and  simulating  the  handing  off  of  observation  responsibilities 
from  site  to  site. 
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4.3  Observation  Types 

In  the  past,  astronomical  measurements  were  limited  to  optical 
sightings  of  celestial  bodies  involving  angular  measurements  and  sighting 
times.  Angular  sighting  observation  pairs  include  azimuth  and  elevation, 
right  ascension  and  declination,  and  meridian  crossing  and  elevation  angle. 
More  recent  astronomical  measurements  involve  electromagnetic  waves 
being  transmitted  through  space.  Electromagnetic  signals  may  be  reflected, 
transponded,  or  transmitted  by  man-made  bodies.  In  addition,  natural  bodies 
are  also  used  to  reflect  radar  signals.  These  more  recent  observation 
techniques  can  provide  precise  interferometric  angular  information  as  well  as 
accurate  range  and  range  rate  measurements.  Independently  of  the 
observation  method  used,  PEP  determines  the  signal  path  from  send  to 
receive  time  using  the  light  time  iterations  described  in  Section  4.1.  PEP 
models  any  effects  which  could  bend  or  distort  this  signal  path  based  upon  the 
type  of  observation.  The  aberration  of  light,  the  Doppler  shift  in  frequency, 
atmospheric  refraction,  ionospheric  distortions,  general  relativity,  and 
interplanetary  plasma  effects  are  all  modeled  in  PEP  to  determine  the  proper 
signal  path.  The  adjusted  signal  path  between  sending  and  receiving  sites,  q, 
is  then  used  to  recreate  observations. 


4.3.1  Azimuth-Elevation  Observations 

Early  astronomical  sightings  measured  the  angles  between  the  apparent 
line-of-sight  to  the  target  body  and  a  reference  frame.  Azimuth  and  elevation 
angle  observations  at  the  observing  site  use  the  vector  normal  to  the 
observing  site  and  the  plane  tangent  to  it  to  describe  the  location  of  the 
observed  body.  The  elevation  angle  is  the  angle  between  the  line-of-sight  and 
its  projection  in  the  tangent  plane.  The  azimuth  angle  is  defined  as  the  angle 
in  the  tangent  plane  between  north  and  the  line-of-sight7  s  projection.  These 
vectors  and  angles  are  depicted  in  Figures  4.3.1-1  and  4.3.1-2. 

In  earth-fixed  Cartesian  coordinates,  the  vector  p  is  [0  0 1]T  and  the  unit 
normal,  n,  is  defined  by  the  longitude  0  and  geodetic  latitude  <(>  as 
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A 

n- 


cos  8  cos  <p 
sin  6  cos  0 


sin0 


(4.3.1-1) 


PEP  transforms  these  vectors  to  the  inertial  coordinate  frame  in  which  q  is 
calculated  using  the  transformations  of  Appendix  C.2. 


Figure  43.1-1:  Azimuth-Elevation  Angle  Vectors  in  die  Meridian  Plane 

The  vector  m,  which  points  along  the  meridian  from  the  observation  site 
towards  the  north,  serves  as  the  north  reference  direction.  This  vector  is 
defined  by  the  following  equation: 

m  =  unit[p  -  (p  •  n)n]  (4.3.1-2) 


Figure  43.1-2:  Azimuth-Elevation  Angles  and  die  Tangent  Plane 


The  projection  of  the  signal  path,  q,  in  the  observation  site's  tangent  plane, 
qp,  is  used  to  define  the  azimuth  angle  and  is  obtained  by  the  equation 
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qp=q-(q*  n)n  (4.3.1-3) 

Additionally,  Figure  4.3. 1-1  shows  the  projection  of  the  signal  path  in  the 
meridian  plane,  qm. 


The  following  equations  are  used  to  calculate  the  observed  azimuth 
and  elevation  angles.  Two-argument  arctangent  routines  should  be  used  for 
improved  accuracy  and  to  distinguish  the  proper  quadrants  [8]. 


elevation  =  sin 


azimuth  =  tan'1 


(4.3. 1-4) 


(4.3.1-5) 


4.3.2  Right  Ascension-Declination  Observations 


When  an  object  is  observed  against  a  star  background,  right  ascension 
and  declination  angles  can  be  determined  from  the  object's  relationship  to 
catalogued  stars.  Right  ascension  and  declination  are  angles  referred  to  an 
inertial  Cartesian  coordinate  system  centered  within  the  observing  body  or 
the  true  equinox  and  equator  of  date.  If  the  line-of-sight  vector  q  is  given  in 
the  relevant  Cartesian  coordinate  system,  [x,  y,  z,]T,  the  right  ascension  and 
declination  angles  are  calculated  from  the  following  relations,  where  once 
again  two-argument  arctangent  routines  should  be  used  for  improved 
accuracy  and  quadrant  determination  [8]. 


declination  =  sin  1 


V 


right  ascension  =  tan"1 


V 

Kx<\  J 


(4.3.2-1) 


(4.3.2-2) 


69 


Figure  4 .3 .2-1:  Right  Ascension  and  Declination  Angies 


4.33  Meridian  Crossing-Elevation  Angle  Observations 

Some  of  the  earliest  observations  recorded  the  time  and  elevation  of  a 
celestial  body's  passage  across  the  earth's  meridian.  These  observations  were 
performed  by  constraining  the  sighting  instrument  in  the  plane  of  the 
meridian  and  measuring  the  elevation  angle  as  the  observed  body  crossed  the 
plane.  Specifying  the  time  of  meridian  crossing  uses  die  rotating  earth  to 
determine  one  of  the  angular  components  of  the  line-of-sight  from  the 
observing  body. 


Figure  4.33-1:  Elevation  Angle  at  Meridian  Crossing 

This  method  is  similar  to  the  azimuth-elevation  observation  with  the 
elevation  angle  constrained  to  the  meridian  plane.  It  is  also  related  to  the 
right  ascension-declination  observation  because  the  local  sidereal  time  at 
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meridian  crossing  converted  to  radians,  is  the  right  ascension.  Based  upon 
the  observing  body's  orbit  and  rotation,  meridian  crossing  time  and  elevation 
angles  are  converted  to  right  ascension  and  declination  angles  [8]. 

4.3.4  Satellite  Look  Angle  Observations 

When  a  satellite  is  the  observing  site,  PEP  defines  the  azimuth  and 
elevation  angles  in  the  satellite's  pitch-roll-yaw  coordinate  system.  The  yaw 
axis,  y,  points  from  the  satellite  to  the  center  of  the  body  which  it  orbits.  The 
roll  axis,  r ,  lies  in  the  orbital  plane  normal  to  the  yaw  axis,  making  an  acute 
angle  with  the  satellite's  velocity  vector.  The  pitch  axis,  p,  is  normal  to  the 
orbital  plane  and  completes  the  right  hand  coordinate  system.  This 
coordinate  frame  is  defined  for  a  lunar  satellite's  position  and  velocity  vectors 
by  the  relations 

y  =  -unit[Xbl],  p  =  «wf[xw x£w],  r  =  yxp  (4.3.4-1) 
~ir 


\ 

\ 

q 

y  I 

Figure  4.3.4-1:  Satellite-Based  Look  Angles 

From  this  coordinate  system,  PEP's  azimuth  angle  is  defined  as  the 
angle  about  the  pitch  axis  and  its  elevation  angle  is  the  rotation  angle  from 
this  point  to  the  line-of-sight.  These  angles  are  depicted  in  Figure  4.3.4-1  and 
the  formulas  for  PEP's  azimuth  and  elevation  angles  are  given  below  [8]. 

(4.3.4-2) 


elevation  =  sin-1 

l 
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azimuth  =  tan-1 


/  _  ~  \ 

iiL 


(43.4-3) 


4.3.5  Interferometry  Observations 

PEP  can  also  determine  angular  measurements  from  interferometer 
measurements.  The  angle  between  the  incoming  wave  and  the  line 
connecting  the  two  observing  sites  can  be  determined  when  two  sites  receive 
an  electromagnetic  signal  and  determine  its  arrival  time  difference. 


If  die  signal  source  (such  as  a  radio  star)  is  sufficiently  far  from  two 
receiving  sites  separated  by  a  distance  d,  then  the  two  signal  paths  are 
considered  parallel.  If  the  signal  arrives  at  site  2  at  a  time  At  after  it  reaches 
site  I,  then  this  signal  has  traveled  the  additional  distance  cAt,  where  c  is  the 
speed  of  light.  The  angle  y  is  then  determined  from  the  relationship 

y^COS'^^j  (43.5-1) 

The  theoretical  value  of  an  interferometer  measurement  of  a  satellite 
is  the  difference  between  the  light  times  from  the  two  observing  sites  at  the 
same  receive  time  and  the  satellite  at  the  signal  transmission  times. 
Detecting  smaller  At's  ».  .^./easing  the  separation  between  the  observing 
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sites  to  even  intercontinental  distances  increases  the  accuracy  of 
interferometric  observations.  The  use  of  hydrogen  maser  atomic  clocks 
allows  angular  measurements  with  milli-arc  second  accuracy.  Previously 
discussed  angular  measurements  are  limited  to  approximately  an  arc  second 
of  accuracy.  PEP  is  coded  to  simulate  and  process  these  observations  and 
model  the  bias  between  the  clocks  at  the  two  receiving  sites  [81. 

4.3.6  Range  Observations 

The  range  between  two  bodies  can  be  determined  from  the  time  delay 
of  a  signal  sent  between  the  bodies.  Range  measurements  depend  upon  a 
precise  knowledge  of  when  the  signal  was  sent  and  when  the  signal  was 
received.  Since  it  is  very  difficult  to  synchronize  two  separated  clocks,  one¬ 
way  ranging  measurements  are  not  often  used,  except  for  multi-GPS  satellite 
observations  where  receiving  site  clock  error  is  measured.  A  two-way  signal 
provides  a  more  accurate  single  satellite  range  measurement  since  a  single 
clock  is  used  to  measure  the  time  delay. 

Two-way  range  observations  can  be  obtained  by  sending  a  pulsed 
electromagnetic  signal  between  two  bodies.  The  time  it  takes  the  pulse  to 
return  to  the  sending  site,  less  any  known  delays,  reveals  the  range  between 
the  bodies  calculated  rigorously  from  light  time  iterations.  Since  the 
electromagnetic  signal  travels  at  the  speed  of  light,  the  range  is  calculated 
from  the  compensated  time  delay.  At,  as 


This  two-way  range  observation  can  also  be  obtained  with  the 
transmission  of  a  continuous  sine  wave  signal  by  shifting  the  phase  of  the 
sine  wave  180°  at  a  certain  repetition  interval  according  to  a  coded  pattern. 
Since  this  phase  shifting  is  preserved  as  the  signal  passes  through  space, 
transponder  electronics,  or  is  reflected  off  of  an  observed  body,  the  receiving 
site  can  recreate  the  phase  shifting  coded  pattern  from  the  received  signal.  By 
correlating  the  sent  coded  pattern  with  the  received  pattern  over  time,  the 
receiving  site  can  determine  the  signal's  round  trip  time  delay.  Using  this 
phase  shift  keying  technique  on  a  continuous  sine  wave  signal  allows  the 
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signal  to  be  used  for  both  range  and  range  rate  observations,  as  well  as  for 
communications,  if  desired  [45]. 


Two-way  time  delay  measurements  in  PEP  are  obtained  directly  from 
knowing  the  site  and  body  positions  at  the  appropriate  times  plus  any 
transponder  delays.  Transponder  time  delays  can  be  determined  prior  to 
launch  or  can  be  estimated  from  the  measurements.  Typically  both  methods 
are  used  by  separating  the  time  delay  into  known  (previously  measured)  and 
unknown  (to  be  estimated)  parts. 


4.3.7  Range  Rate  Observations 

Range  rate  observations  measure  the  rate  at  which  the  observation  site 
and  observed  body  are  approaching  one  another.  These  observations  take 
advantage  of  the  Doppler  effect  upon  electromagnetic  signals.  Since  this 
observation  measures  the  change  in  frequency  between  send  and  receive 
time,  it  has  the  same  difficulty  with  one-way  observations  as  range 
measurements.  Coherent  two-way  signal  paths,  however,  provide  very 
accurate  measurements  of  the  range  rate  between  two  bodies. 

The  Doppler  effect  is  a  shift  in  the  frequency  of  an  electromagnetic 
wave  radiated,  reflected  or  received  by  an  object  in  motion.  This  frequency 
shift  is  a  result  of  the  expansion  or  compression  of  electromagnetic  waves 
along  the  direction  of  a  moving  source  [43].  This  compression  in  the 
direction  of  motion  is  visualized  in  Figure  4.3.7-1  below.  For  astronomical 
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range  rate  measurements,  the  frequency  of  the  signal  sent  from  site  1  to  site  2 
is  shifted  due  to  the  relative  velocity  of  site  1  at  send  time  and  site  2  at  receive 
time  along  the  line-of-sight.  From  the  light  time  iteration  process  mentioned 
in  Section  4.1,  PEP  determines  this  relative  velocity  from  the  state  of  the  two 
sites  at  the  two  times  in  an  inertial  reference  frame. 


Figure  4. 3. 7-1:  Waves  Radiated  from  a  Stationary  and  Moving  Source 

For  a  range  rate  of  R,  the  one-way  Doppler  shift  of  a  frequency  source, 
fs,  is  given  by  the  formula 

Af  =  -fs~  (43.7-1) 

R  is  positive  when  the  bodies  are  moving  apart  and  the  signal's  frequency  is 
decreasing  due  to  the  Doppler  shift.  This  formula  is  used  to  determine  the 
frequency  received  or  reflected  at  the  observed  body.  For  radio  wave 
transponder  observations,  the  received  frequency  is  multiplied  by  a  rational 
factor  q,  so  that  the  body  transponds  a  different  frequency,  f$2,  from  the  signal 
it  receives. 


/.2=<l/S 


( 


1-5 

c 


(43.7-2) 


Once  again,  after  transmission,  this  frequency  shifts  due  to  the  Doppler  effect 
by  (4.3.7-1)  where  fs  is  the  send  frequency  from  (43.7-2)  and  the  range  rate 
depends  upon  the  new  send  and  receive  time  positions  and  velocities. 

For  approximately  instantaneous  send  and  receive  times  with  no 
transponder  frequency  shifting  (q  =  1),  the  round  trip  Doppler  shift  is 
approximated  by  the  single  expression  [43] 
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D 

A/  =  -2/s  —  (4.37-3) 

The  NASA  Deep  Space  Network  (DSN)  measures  a  coherent  count  of 
the  zero  crossings  of  the  Doppler  signal/  refered  to  as  Integrated  Doppler.  A 
DSN  site  sends  a  continuous  sine  wave  signal  towards  a  spacecraft  at  a 
frequency  controlled  by  an  atomic  frequency  standard.  The  receiving  DSN 
station  can  be  the  same  or  different  from  the  sending  site  because  of  the 
network's  precise  synchronization.  The  receiving  station  pre-multiplies  the 
transmitted  frequency  by  the  spacecraft's  transponder  translation  factor  and 
then  subtracts  this  frequency  from  the  received  signal.  Any  system  biases  are 
then  added  to  this  differenced  frequency/  if  known/  or  they  can  be  estimated 
in  the  fitting  process.  A  counter  at  the  receiving  site  is  incremented  for  every 
positive  traveling  zero  crossing  of  the  differenced  frequency.  This  counter  is 
read  at  uniform  intervals/  At,  to  determine  the  interval's  Doppler  frequency 
shift  in  cycles  per  second  (Hz).  Additionally,  a  resolver  is  used  to  give 
fractional  cycle  resolution  of  the  zero  crossings  [36]. 

Within  PEP,  the  theoretical  value  of  this  observable  is  the  difference  in 
round  trip  phase  delays  at  the  count  starting  and  ending  receive  times 
multiplied  by  the  sending  frequency.  PEP  also  accounts  for  any  transponder 
frequency  translation  when  calculating  the  theoretical  value  of  the 
observation.  The  exact  formula  for  the  Doppler  count  observable  is  coded  in 
PEP  [8]  and  is  approximately  equal  to  the  instantaneous  Doppler  shift  over  the 
count  interval. 


4.4  Partial  Derivatives  of  Satellite  Observations 

The  partial  derivatives  of  the  theoretical  value  of  an  observation  are 
required,  along  with  the  observation  residuals,  to  estimate  orbital  initial 
conditions  and  parameters.  The  theoretical  value  of  an  observation  is  a 
function  of  the  receiving  site  coordinates  at  signal  receive  time,  the  observed 
body's  coordinates  at  transpond  time,  the  sending  site's  coordinates  at  send 
time,  and  other  parameters.  If  bent  pipe  observations  are  used,  the  theoretical 
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value  is  also  a  function  of  any  additional  body's  coordinates  at  the  appropriate 
times.  For  a  given  observation,  T,  this  functional  relationship  is  denoted  by 

r<fr) = /(x;„  (tr  >,x;t  a,  ),xi  (t, ) . P)  (4.4-d 

where  t,  is  the  receive  time,  r.s.  is  the  receiving  site,  t.s.  is  the  transponding  or 
reflection  site,  t*  is  the  transpond  or  reflection  time,  s.s.  is  the  sending  site,  and 
t,  is  the  send  time.  The  ellipses  denote  that  T  may  be  a  function  of  other 
coordinates  at  other  times,  depending  upon  the  observation  method,  ft  is  the 
vector  of  parameters  besides  motion  which  affect  r. 

As  the  theoretical  observations  are  determined,  their  partial 
derivatives  with  respect  to  the  coordinates  of  interest  can  also  be  calculated. 
Depending  upon  the  type  of  observation,  PEP  also  calculates  the  partial 
derivatives  of  the  observation  with  respect  to  the  parameters,  /?,  which  affect 
the  observation,  such  as  measurement  biases  and  transponder  delays.  The 
partial  derivative  of  the  theoretical  observation,  T,  with  respect  to  a  parameter 
of  interest,  a,  is  then  calculated  by  the  chain  rule. 

dr  dr  <?x;,  dr  <?x*  dr  ax',  ar  dp ,  „ 

da  dx;s  da  dx;3  da  dX^  da  dp  da 

The  parameter  a  would  be  any  parameter  to  be  estimated,  such  as  orbital 
initial  conditions,  gravitational  harmonic  coefficients,  observing  site 
coordinates  or  observation  biases.  The  partial  derivative  of  a  satellite's 
coordinates  are  determined  by  interpolation  from  the  satellite's  integration 
file  using  the  methods  discussed  in  Section  3.1.  The  partial  derivative  of  a 
site  coordinate  is  a  function  of  the  site's  spherical  or  cylindrical  coordinates. 

Once  calculated,  these  partial  derivatives  are  written  to  an  observation 
output  file  for  each  receive  time  for  each  observation  type  in  an  observation 
series.  The  theoretical  observations,  their  partial  derivatives,  and  the 
observation  residuals  for  each  observation  series  on  the  output  file  can  then 
be  used  to  calculate  adjustments  to  the  parameters  a. 
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Chapter  Five 


Parameter  Estimation 
and 

Prediction  of  Uncertainty 


5.1  Parameter  Estimation 

Estimation  is  the  process  of  extracting  information  from  a  collection  of 
observations  (measurements)  to  develop  a  better  understanding  of  the 
observed  behavior.  For  this  thesis,  the  parameter  estimation  process  attempts 
to  obtain  the  best  set  of  lunar  gravitational  harmonic  coefficients.  The  values 
obtained  in  the  estimation  process  are  best  in  the  sense  that  the  estimated 
model's  generated  observations  provide  the  closest  match  with  the  "truth" 
model  simulated  observations.  This  does  not  mean  that  a  better  fit  could  not 
be  realized  with  more  or  different  "truth"  model  observations.  It  also  does 
not  imply  that  a  better  fit  could  not  be  obtained  with  a  better  model  of  the 
behavior.  The  estimation  process  additionally  provides  insight  into  the 
uncertainty  of  the  estimates'  ability  to  model  the  behavior.  This  information 
can  be  used  to  analyze  and  evaluate  the  estimated  model. 


5.1.1  Observation  Vector  as  a  Random  Variable 

Using  the  gravitational  "truth"  model,  observations  were  generated  for 
several  proposed  sensing  schemes.  These  observations  are  dependent  on  the 
lunar  gravitational  field  and  orbital  states  of  the  satellites.  Together  these 
observations  form  the  N-dimensional  vector  z  where  N  is  the  total  number 
of  measurements. 
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With  initial  guesses  for  1)  the  lunar  gravitational  field  coefficients, 
2)  satellite  initial  conditions  (initial  osculating  orbital  elements  for  each 
satellite  used  in  the  sensing  scheme),  and  3)  measurement  biases,  theoretical 
values  of  the  estimated  model  observations,  the  difference  between  the  two 
sets  of  observations,  and  the  partial  derivatives  of  the  estimated  model 
observations  are  generated.  The  estimated  model  observations'  theoretical 
values  are  a  function  of  the  observation  time,  the  estimated  model 
parameters,  and  other  factors,  such  as  the  positions  of  the  observing  sites. 
These  other  factors  are  treated  as  known  quantities  for  the  estimation  process. 
Let  ag  be  the  vector  of  unknown  gravitational  parameters  and  ft  be  the 
vector  of  unknown  measurement  biases  and  any  other  unknown  parameters. 
For  a  sensing  scheme  with  N&  satellites,  there  are  Nj,  vectors  an  of  unknown 
orbital  initial  conditions.  If  a  is  the  vector  of  all  of  these  unknown  system 
parameters,  then  the  theoretical  observation  at  time  t  is  a  function  of  t  and  a. 


si=[h 
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(5.1.1-3) 

Due  to  imperfect  knowledge  of  the  gravitational  field  coefficients  and 
satellite  initial  conditions  and  the  difference  between  the  fit  and  truth 
models,  the  theoretical  observations  do  not  perfectly  recreate  the  "truth" 
model  observations.  For  each  measurement  time,  t,  there  is  a  measurement 
error,  0(t),  between  the  "truth"  model  and  theoretical  observations.  With  the 
"true"  values  for  the  parameters  a,  the  remaining  measurement  error,  0,  is 
assumed  to  be  a  zero-mean  random  variable  (disregarding  the  difference 
between  the  fit  and  truth  models).  Since  6  is  the  result  of  several 
independent  random  causes,  it  will  be  assumed  to  be  normally  distributed  by 
the  central  limit  theorem.  The  collection  of  N  measurement  errors  over  the 
observation  period  is  then  6. 

From  the  previous  definitions,  it  follows  that  the  vector  of 
measurements  is  a  random  vector  composed  of  a  deterministic  value,  the 
theoretical  observation,  and  a  random  quantity,  the  measurement  error. 
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1 

or  (5.1. 1-4) 

z  =  f(t-,a)  +  6 

Since  the  measurement  error,  0,  is  assumed  normally  distributed,  the 
measurement  vector  is  also  normally  distributed  with  the  following  statistics: 

fi  =  E{z}  =  E[f(t;a)}  +  E{9}  =  f(t;a)  (5.1.1-5) 

E{(z  -  fi)(z -fiY}  =  E{00T}  =  ©  (5.1.1-6) 


5.1.2  Maximum  Likelihood  Estimation 

The  maximum  likelihood  estimate  selects  the  parameters  a  such  that 
the  probability  of  z,  the  likelihood  function,  is  maximized.  This  method 
estimates  parameter  values  so  that  the  "truth"  model  observations  are  the 
most  likely  ones  to  have  occurred.  Since  z  is  a  normal  random  variable,  its 
probability  density  is  the  normal  joint  probability  density: 

p(z;«)  =  /(2  J)w|e|expH(g  -/(*;“))]  e-i-2-i) 

The  maximum  likelihood  estimate  maximizes  the  likelihood  function 
(5.1.2-1)  and  minimizes  the  negative  log-likelihood  function  below: 

£(z;a)  =  -ln[p(z;a)]  '5.1.2-2) 

£(z;  a)  =  ln(2«)  +  |ln(|0|)]  +  ±(z-  f(t;  a)f  0'1  (z  -  /(f;  a)) 

(5.1.2-3) 

The  maximum  likelihood  estimate  a  therefore  satisfies  the  equation 
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(5. 1.2-4) 


The  maximum  likelihood  estimate  &  is  obtained  starting  with  an 
initial  guess  of  the  parameters,  a0.  Based  upon  the  theoretical  observations 
generated  with  a0,  adjustments.  Ad,  to  the  guessed  parameters  are  sought 
such  that 

A  a  =  a-a0  (5. 1.2-5) 


Expanding  Equation  (5.1.2-4)  in  a  Taylor  series  expansion  about  a0  yields 


Aa  +  0(Aa2)  (5.12-6) 


The  parameter  adjustments,  A  a,  are  determined  by  taking  the  partial 
derivatives  of  the  negative  log-likelihood  function  from  Equation  (5. 1.2-3) 
(assuming  that  the  term  in  brackets  does  not  depend  on  the  parameters  a) 
and  by  neglecting  the  higher  order  terms  in  die  Taylor  series  expansion.  The 
second  derivative  of  the  negative  log-likelihood  function  with  respect  to  a  is 
replaced  with  its  expected  value,  which  can  be  theoretically  expressed  as  the 
expected  value  of  the  dyadic  product  of  first  partial  derivatives  [42]: 

E  W  y  i  =  J(  y 

day  da  )  y  da  J  da 


Assuming  that  0  does  not  depend  on  a,  Equations  (5.1.2-3)  and  (5.1.2-7)  imply 


n\±{ 

da  J 

«.  - 


(5.1. 2-8) 


Replacing  the  Hessian  of  second  partial  derivatives  in  Equation  (5.1.2-6)  by  its 
expected  value  (the  Fisher  information  approximation)  leads  to  the  linear 
matrix  equation,  called  the  normal  equations,  for  the  adjustments  to  the 
parameters.  These  equations  are  formed  with  the  partial  derivatives  of  the 
theoretical  observations  with  respect  to  the  estimated  parameters,  the 
observation  residuals,  and  the  measurement  error  statistics. 
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AAa  =  b 

or 


r  df(z;a0))'J  <if(z;a0) 
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(5.1.2-9) 


5.13  Weighted  Least  Squares  Estimation 

Least  squares  estimation  selects  the  parameters  St  so  as  to  minimize  the 
sum  of  the  squares  of  the  deviations  between  the  "truth"  model  and 
theoretical  measurements.  For  N  observations,  the  least  squares  estimate 
seeks  to  minimize  the  quantity 

N 

Q=X(zi -/■(*(;«))  (5.13-1) 

i=l 

If  the  measurements  are  of  varying  quantities  and  units,  and  some 
measurements  are  more  reliable  than  others,  the  weighted  cost  function  Q  is 
used. 


<3=1 

i=i 


(z,  -/,(*■;  «))2 

Wi 


(5. 1.3-2) 


where  W,  is  a  sequence  of  weighting  values  [21].  A  sequence  of  observations 
may  involve  a  wide  variety  of  physical  quantities,  i.e.  distances,  angles, 
temperatures,  frequencies.  Typically  the  error  weightings  non- 
dimensionalize  these  disparate  measurements  in  the  cost  function. 

The  least  squares  or  weighted  least  squares  estimate  is  found  by  setting 
the  derivative  of  the  cost  function  with  respect  to  the  estimation  parameters 
a  equal  to  zero. 
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The  weighted  least  squares  and  maximum  likelihood  methods  provide 
the  same  estimates  if  the  measurement  residuals,  9,  are  uncorrelated. 


e{<M>,}=o 


V  i*j 


(5.1.4-1) 


By  assuming  that  this  measurement  noise  is  uncorrelated,  the  covariance 
matrix,  6  defined  in  Equation  (5.1.1-6)  becomes  the  diagonal  matrix 


of  0  •  •  0 
0  of  •••  0 


(5. 1.4-2) 


L  0  0  o-NJ 

where  a\  is  the  variance  and  an  is  the  standard  deviation  of  the  nth 
measurement.  Inserting  this  diagonal  matrix  into  Equation  (5. 1.2-1)  and 
(5.1.2-3)  for  the  likelihood  and  negative  log-likelihood  functions  results  in  the 
following: 


p(z;a)  =  — g - exp 

(27t)2  <V”0-n 


,  y(z, 


(5. 1.4-3) 


C(z; a )  =  [f ln(2;r)  +  ln^ •  •  •  crN )]  ^ +  (5>1.4^) 


Both  the  maximum  likelihood  and  weighted  least  squares  estimates 
are  determined  by  setting  the  partial  derivative  of  a  cost  function  with  respect 
to  a  equal  to  zero.  Since  the  maximum  likelihood  estimate  assumes  that  the 
constant  part  in  brackets  [  1  from  (5.1.4-4)  does  not  depend  upon  the 
parameters  a,  the  two  methods  (compare  Equations  (5.1.3-2)  and  (5.1.4-4)) 
provide  identical  estimates  when  uncorrelated  measurement  error  variances, 
o\,  are  used  to  weight  the  measurement  residuals.  Since  PEP's  parameter 
estimation  routine  assumes  uncorrelated  measurement  errors,  its  estimation 
technique  is  referred  to  as  least  squares  maximum  likelihood  estimation. 
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5.1.5  Incorporation  of  A  Priori  Information 

Frequently,  existing  models  provide  a  priori  information  about  some 
of  the  estimation  parameters.  The  estimation  process  can  be  shortened  or 
simplified  by  including  this  a  priori  information  in  the  normal  equations. 

Suppose  a  priori  information  has  estimated  the  value  of  the  parameter 
a,  as  a{  with  an  uncertainty,  or  standard  deviation,  of  d,.  These  a  priori 
parameter  values  are  then  grouped  into  the  m-dimensional  vector  a,  where 
m  is  the  total  number  of  parameters  estimated.  A  zero  value  is  assigned  to 
any  parameters  when  no  a  priori  knowledge  is  available.  The  variances  of 
the  a  priori  estimates  are  collected  into  the  m  x  m  diagonal  matrix  £.  For  any 
parameters  without  a  priori  standard  deviation  information,  their  diagonal 
element  is  infinity,  although  in  practice  some  reasonably  large  value  will 
suffice.  A  full  covariance  matrix  can  be  used  for  £  if  correlations  are 
available  for  the  a  priori  parameter  estimates  a. 

With  a  priori  information,  instead  of  minimizing  Equation  (5.1. 2-3), 
the  following  functional  is  minimized  [11]. 

Q  =  (z-/(f;  a))1©'1  (z  -  f(t ;a))  +  (a  -  a)T  XT1  (a  -  a)  (5.1.5-1) 

For  parameters  where  no  a  priori  standard  deviation  information  is  available, 
the  cost  function  remains  the  same  as  in  (5. 1.2-3). 


Setting  the  partial  derivative  of  this  cost  function  with  respect  to  the 
parameters  a  equal  to  zero  and  making  the  same  assumptions  as  in  Section 
5.1.2,  yields  the  following  linear  matrix  equations. 


A  = 


df(z;a0)  . 
da 


(5. 1.5-2) 


b  = 


r  df(z}g0)^ 
da 


0"1  (z  -  f(z ;  a0 ))  -  XT1  (a0  ~  «)  (5- 1.5-3) 


such  that 

A  a  =  A-1  b 


(5. 1.5-4) 
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From  studies  of  the  earth's  gravitational  held,  William  Kaula  proposed 
an  empirical  relationship  for  the  variance  of  the  gravitational  harmonic 
coefficients  [24].  Estimates  of  the  lunar  gravitational  harmonic  coefficients  led 
him  to  believe  that  the  same  relationship  was  valid.  For  the  lunar 
coefficients,  Kaula  adjusted  the  constant  coefficient  to  account  for  the 
gravitational  field  strength  difference  between  the  moon  and  earth  [37].  From 
their  studies  of  the  lunar  gravitational  field.  Bills  and  Ferrari  chose  a  slightly 
different  semi-empirical  formula  for  the  covariance;  one  which  had  been 
proposed  by  Vening-Meinesz.  After  having  difficulty  converging  upon 
gravitational  harmonic  coefficients.  Bills  and  Ferrari  included  the  Vening- 
Meinesz  a  priori  covariance  information  and  developed  a  16  x  16  lunar 
spherical  harmonic  potential  model  with  coefficient  uncertainties  [12]. 
Research  scientists  at  JPL  continue  to  use  this  method  to  develop  50  X  50,  and 
even  75  X  75  spherical  potential  models  of  the  moon  [27]. 

The  method  of  incorporating  a  priori  information  into  the  normal 
equations  was  not  used  for  these  simulations.  This  thesis  evaluates 
gravitational  sensing  schemes  which  may  be  employed  to  develop  gravity 
field  models.  As  such,  an  arbitrary  lunar  gravitational  "truth''  model  was 
developed  and  there  was  no  a  priori  information  available  regarding  this 
"truth"  model.  When  a  gravitational  sensing  method  is  selected  and  its 
mission  flown,  a  priori  lunar  gravitational  field  information  can  and  should 
be  used  in  the  development  of  new  gravitational  field  models. 


5.1.6  Solution  to  the  Normal  Equations 

Since  the  normal  equations  [(5.1.2-9)  or  (5.1.5-4)]  are  linear,  the 
adjustments  to  the  parameters  can  be  solved  by  various  numerical 
techniques.  PEP  uses  the  Gauss-Jordan  method  to  simultaneously  solve  the 
normal  equations  and  invert  the  coefficient  matrix.  This  method  uses 
diagonal  pivots  without  interchanging  rows  or  columns,  so  that  only  the 
lower  diagonal  half  of  the  symmetric  coefficient  matrix  A  is  stored  in 
memory. 

Numerical  problems  can  arise  if  the  coefficient  matrix  is  ill- 
conditioned.  This  could  happen  if  the  vector  of  parameter  adjustments  A  a 
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consisted  of  widely  different  sized  quantities  or  quantities  with  widely 
different  units.  To  prevent  these  numerical  difficulties,  PEP  scales  the 
normal  equations  prior  to  solving  them.  This  process  scales  each  row  and 
column  of  the  symmetric  coefficient  matrix  by  the  square  root  of  the  diagonal 
element.  To  preserve  the  equality,  each  row  of  the  right  hand  side  vector  is 
also  scaled  by  the  same  factor.  After  the  adjustments  to  the  parameters  sure 
solved,  PEP  unscales  its  rows  to  provide  the  proper  units  and  values  to  the 
parameter  adjustments  [8]. 

This  scaling  process  can  be  avoided  by  using  the  square  root  of  the 
normal  equations  [11].  This  method  takes  advantage  of  the  properties  of  the 
symmetric  coefficient  matrix  and,  by  using  the  square  root  of  the  normal 
equations,  lessens  the  effect  of  disparate  units  or  scale  factors.  The  JPL  orbit 
fitting  software  uses  this  method  [27]. 

Since  the  equation  for  Ad  was  obtained  by  neglecting  second  and 
higher  order  terms  in  a  Taylor  series  expansion,  this  adjustment  will  not 

A 

yield  a  exactly.  An  iterative  technique  must  be  used  to  approach  a  maximum 
likelihood  estimate  for  the  parameters.  Once  the  parameter  adjustments  are 
determined,  the  initial  parameter  guesses  are  adjusted.  The  equations  of 
satellite  motion  are  then  re-integrated  and  the  theoretical  observations, 
partial  derivatives,  and  residuals  are  recalculated.  The  normal  equations  are 
then  reformed  and  a  second  set  of  parameter  adjustments  are  determined. 
These  iterations  continue  until  the  process  converges. 


5.1.7  Other  Estimation  Techniques 

For  some  applications,  least  squares  maximum  likelihood  may  not  be 
the  optimum  method  for  estimating  the  gravitational  coefficients.  If  noise 
were  included  in  the  satellite  dynamics,  Kalman  filtering  techniques  would 
be  more  appropriate.  A  linearized  Kalman  filter  about  a  nominal  satellite 
trajectory  or  an  extended  Kalman  filter  for  which  the  reference  trajectory  is 
updated  for  each  observation  could  be  used  [31,  Chapter  9].  The  Kalman  filter 
would  then  estimate  the  gravitational  parameters  as  well  as  the  satellite's 
state  by  augmenting  the  state  vector  to  include  them  as  non-dynamic  state 
variables. 
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Since  Kalman  filtering  techniques  become  extremely  cumbersome  as 
higher  degree  models  are  estimated,  a  more  efficient  technique  may  be 
necessary.  Maximum  likelihood  system  identification  combines  the  two 
approaches  by  performing  a  Kalman  filter  on  the  satellite  state,  X*,  and  a 
maximum  likelihood  estimate  on  the  satellite's  orbital  initial  conditions  and 
the  gravitational  harmonic  coefficients  [31,  Chapter  10].  In  this  method,  the 
satellite's  motion  is  propagated  with  a  Kalman  filter  and  then  a  maximum 
likelihood  adjustment  is  made  to  the  initial  conditions  and  parameters.  The 
process  is  then  repeated  until  convergence. 

Unmodeled  forces  such  as  radiation  pressure,  fuel  tank  leakage,  and 
higher  order  gravitational  harmonic  effects  will  cause  noise  in  the  satellite 
dynamics.  Because  of  this  noise,  one  of  these  techniques  should  be  used. 
Least  squares  maximum  likelihood  estimation,  however,  was  sufficient  for 
this  thesis'  evaluation  of  sensing  methods. 


5.2  Statistical  Prediction  of  Uncertainty 

Solving  the  normal  equations  produces  information  about  the 
uncertainty  of  the  estimated  parameters  a  as  well  as  the  values  of  the 
parameters  themselves.  The  parameter  estimation  covariance  matrix 
provides  a  measure  of  the  uncertainty  of  the  estimates.  This  additional 
information  can  be  used  to  evaluate  the  resulting  estimated  model  of  the 
observed  behavior. 

The  coefficient  matrix  A  from  (5.1.2-9)  and  (5.1.5-4)  is  the  Fisher 
information  matrix  and  by  the  Cramer-Rao  lower  bound,  the  covariance  of 
any  unbiased  estimate  is  greater  than  or  equal  to  its  element  in  the  inverse  of 
the  Fisher  information  matrix  [42,  46].  Additionally,  maximum  likelihood 
estimates  have  the  desirable  property  that  they  are  asymptotically  unbiased 
with  the  Cramer-Rao  lower  bound  obtained  as  the  number  of  observations 
approaches  infinity  [42].  This  bound  only  applies  if  the  nonlinear  system  has 
been  modeled  correctly  in  (5.1.1-4). 


88 


If  the  assumed  standard  deviation  of  the  measurement  errors,  an<  are 
correct  and  the  modeling  is  correct  for  both  the  equations  of  motion  and  the 
observations,  the  root  mean  square  (rms)  of  the  post-fit  observation  residuals 
divided  by  the  assumed  measurement  errors  would  be  approximately  one: 

nns  =  -  f(t;  a))T0-'(z  -  /(f;a>)  (5.2-1) 

Typically  the  number  of  observations  is  much  greater  than  the  number 
of  parameters,  so  the  l/(N-m)  term  is  replaced  by  1/N.  When  the  above  rms 
differs  greatly  from  one,  the  parameter  covariance  matrix  produced  by  the 
inverse  of  the  coefficient  matrix  from  the  normal  equations  (A-1),  by  the 
Cramer-Rao  lower  bound,  should  be  multiplied  by  this  rms  to  obtain  a  truer 
estimate  of  the  uncertainty  in  the  parameter  estimates.  This  adjustment 
accounts  for  incorrect  values  of  the  measurement  standard  deviations,  a„  in 
8,  but  does  not  account  for  any  modeling  error  effects  on  the  uncertainty. 

The  parameter  covariance  matrix,  Z,  is  then  obtained  as 

I  =  rms  x  A*1 

The  standard  deviation  of  parameter  estimate  a,  is 

Gi  = 

and  the  correlation  between  parameter  estimates  a,  and  &j  is 


PEP  saves  the  matrix  A'1  and  the  rms  of  the  observation  residuals 
divided  by  the  assumed  measurement  errors  resulting  from  an  estimation 
run  so  they  can  be  used  to  predict  the  uncertainty  of  an  orbit  propagated  with 

A 

the  estimated  parameters  a. 

Let  X*  be  the  state  of  a  lunar  satellite  at  time  to  assumed  to  be  known 
perfectly,  perhaps  by  determination  from  navigation  satellite  observations. 
These  satellite  initial  conditions  are  numerically  integrated  in  time  using  the 
estimated  gravitational  model.  Partial  derivatives  of  the  satellite's  motion 
with  respect  to  the  parameters  a  are  also  numerically  integrated.  The  state 


(5.2-2) 

(5.2-3) 

(5.2-4) 
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covariance  matrix,  E,  for  the  uncertainty  in  the  satellite  state,  can  then  be 
calculated  at  any  time  t  as 


da  da 


(5.2-3) 


The  standard  deviation  of  the  uncertainty  in  any  one  coordinate 
direction  is  the  square  root  of  that  diagonal  element.  The  root  sum  square 
(rss)  of  the  position  and  velocity  uncertainties  are  then 

tss(R)  =  -y/Eu  +5^  +  E33  (5.2-4) 

rss(V)  =  y'E^  +S55  +Zm  (5.2-5) 


Since  the  state  covariance  matrix  E  is  based  upon  the  augmented  state 
vector  (3.3-1),  it  is  partitioned  as  follows: 


5^(0  ]  ERV(t) 

ZRV(i)tz™(i) 


(52-6) 


Often,  the  cross  correlations  between  position  and  velocity  ERV  are  neglected 
and  the  covariance  matrix  for  position  ERR  and  velocity  Evv  are  used 
separately  to  evaluate  the  uncertainty  of  orbit  prediction.  The  state 
uncertainty  analyses  performed  for  the  estimated  lunar  gravity  fields  in 
Chapter  Seven  use  these  individual  state  covariance  matrices  rather  than  the 
augmented  state  covariance  matrix. 

The  inertial  (x,y,z)  coordinate  uncertainties  thus  obtained  can  also  be 
converted  to  a  local  vertical,  local  horizontal  coordinate  frame  to  provide 
navigation  uncertainty  information.  The  local  vertical,  local  horizontal 
frame  is  defined  by  the  vertical  (VT),  down  range  (DR),  and  cross  track  (CT) 
directions.  Their  translations  from  inertial  coordinates  are  obtained  by 


uVT(t)  =  unit[X(t )] 

(5.2-7) 

Mcj  (o  =  unit  J&* )  x  X(/)] 

(5.2-8) 

^DR  ^  =  ^VT  W  X  Uqj  ( t ) 

(5.2-9) 
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The  transformation  matrix  from  local  vertical,  local  horizontal  to 
inertial  coordinates  is  then 

Rn.(0  =  ^ ct W  1  (5.2-10) 

and  this  transformation  can  be  used  to  obtain  the  covariance  matrix  for  the 
uncertainty  in  the  satellite  state  in  the  LVLH  frame. 

=  Rj. <oSf*  ( t )  Ril(0  (5.2-11) 

E  lv(t)  =  R*  0)5™  ( t )  R^U)  (5.2-12) 

The  standard  deviation  of  the  uncertainty  in  any  one  LVLH  coordinate 
direction  is  the  square  root  of  that  diagonal  element.  With  these  formulas, 
the  uncertainty  of  position  and  velocity  in  the  vertical,  cross  track,  and  range 
directions  can  be  calculated  at  each  time  t  due  to  uncertainties  in  the 
estimated  parameter  Sc.  These  routines  were  coded  in  analytical  software 
separate  from  PEP  and  each  estimated  gravity  model  was  evaluated  in  this 
manner.  This  analysis  was  performed  for  both  the  PEP  supplied  covariance 
matrix  (A*1)  and  the  covariance  adjusted  by  the  rms  of  the  observation 
residuals  divided  by  the  assumed  measurement  errors.  These  analyses 
simulated  the  propagation  of  uncertainties  onto  the  far  side  of  the  moon  after 
obtaining  a  navigation  fix  on  the  near  side. 

The  previous  statistical  uncertainty  prediction  was  performed  using 
the  assumption  that  the  lunar  satellite  initial  state,  X*,  was  known  perfectly. 
If  this  is  not  the  case,  the  same  process  for  calculating  the  uncertainty  of  orbit 
prediction  would  still  apply.  In  general,  the  vector  of  estimated  parameters  Sc 
would  be  augmented  to  include  estimates  of  the  satellite  initial  conditions 
and  the  process  repeated.  Since  this  thesis  is  concerned  with  the  uncertainty 
due  to  a  mismodeled  lunar  gravitational  field,  perfect  knowledge  of  the 
satellite's  initial  conditions  was  assumed. 
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5.3  Prediction  Uncertainty  Due  to  Mismodeling 

The  statistical  uncertainty  analysis  presented  in  the  previous  section  is 
based  upon  the  assumption  that  the  estimated  gravitational  model  correctly 
models  the  lunar  gravity  field.  Due  to  the  inclusion  of  mascons  in  the 
"truth"  model,  the  finite  spherical  harmonic  estimation  model  could  not 
correctly  model  the  observed  behavior.  In  most  real  world  cases  there  is  no 
way  to  account  for  or  analyze  the  errors  between  modeled  and  true  behavior 
because  the  modeling  errors  cannot  be  separated  from  the  other  errors.  For 
this  thesis,  however,  a  direct  comparison  between  the  true  and  predicted 
behavior  can  be  made,  since  both  the  "truth"  and  fit  models  are  available. 

Given  satellite  initial  conditions,  X*,  orbits  for  the  truth  and  fit  models 
can  be  numerically  integrated.  At  each  point  in  time,  the  position  and 
velocity  deviations  between  the  two  models  can  be  calculated.  These 
modeling  errors  in  position  and  velocity  can  then  be  transformed  from 
inertial  coordinates  to  the  local  vertical,  local  horizontal  frame  using  the 
transformations  (5.2-7),  (5.2-8),  and  (5.2-9).  The  root  sum  square  of  the 
position  and  velocity  errors  due  to  gravitational  field  mismodeling  can  also 
be  calculated.  The  impact  mismodeling  has  upon  navigation  uncertainty 
prediction  is  revealed  by  comparing  the  modeling  errors  with  the  statistical 
uncertainty  predictions. 

Once  again,  using  software  separate  from  PEP,  this  analysis  was 
performed  for  orbits  propagated  from  the  lunar  near  side  to  the  far  side  for 
one  orbital  period.  Additionally,  the  uncertainty  due  to  mismodeling  was 
analyzed  by  comparing  the  "true"  path  of  a  lunar  landing  when  the 
maneuver  was  calculated  based  upon  simulations  using  the  estimated  model 
and  then  executed  in  the  "true"  lunar  gravitational  field.  These  landing 
errors  provide  a  feel  for  the  impact  lunar  gravitational  field  mismodeling 
will  have  upon  the  execution  of  future  space  missions. 
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Orbit  Fitting 
and 

Gravity  Estimation 


6.1  Methodology 

Several  different  lunar  gravitational  sensing  methods  were  evaluated 
using  the  PEP  orbit  integration,  observation  generation,  and  parameter 
estimation  capabilities  described  in  Chapters  Three,  Four  and  Five.  For  each 
sensing  scheme,  the  following  procedures,  depicted  in  Figure  6.1-1,  were  used 
to  generate  "truth"  model  observations  and  estimate  gravitational  field 
coefficients. 

First,  each  sensing  scheme  was  analyzed  to  determine  the  nature  of  the 
satellite  orbits  employed,  the  type  of  observations  generated,  and  the 
observation  accuracy  and  frequency.  Using  the  "truth"  model  and  these 
descriptions  of  the  observation  method,  PEP  simulated  the  mission's  orbits 
and  observations,  including  lunar  occultations.  Auxiliary  software  routines 
converted  the  integration  output  to  selenographic  osculating  orbital  elements 
(a,e,i,fl,ci),M)  as  well  as  selenographic  spherical  coordinates  (r,0,<(>)  versus  time. 
This  data  was  plotted  to  analyze  the  orbit's  stability  over  the  integration 
period  (14  to  28  days).  For  stable  orbits,  the  "truth"  model  observation  file  was 
then  used  to  estimate  the  fit  model  lunar  spherical  harmonic  coefficients  and 
satellite  initial  conditions. 

For  the  estimation  process,  the  gravitational  harmonic  coefficients  and 
satellite  initial  conditions  were  perturbed  from  the  "true"  values  used  to 
generate  the  observations.  The  gravity  field  was  altered  to  reflect  the  fit 
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model's  degree  and  order  ("truth"  model  usually  included  mascons).  From 
the  initial  perturbed  parameters,  a0,  the  satellite  orbits  were  propagated, 
theoretical  values  of  the  observations  were  calculated,  and  the  difference 
between  the  "truth"  model  observations  and  these  theoretical  observation 
values  (observation  residuals)  were  calculated.  Partial  derivatives  of  the 
motion  and  partial  derivatives  of  the  theoretical  observations  were  calculated 
with  respect  to  all  parameters  to  be  estimated.  From  the  observation  residuals 


Truth  Model  & 
Satellite  l.C.‘s 
Input 


PEP  Orbit  I 
Integration! 


Truth  Orbir 
Integration 
.  File  J 


PEP 

Observation 

Generation 


’Truth"  Model  Observations 


Fit  Model  & 
Satellite  l.Cs 
Input 


'Auxiliary 

Analysis 

^Software, 


PEP  Orbit 

r  Orbit  ^ 

Integration 

A 

^  1  integration 
^w/  Partials  FLUj 

Input  for  New 
Parameter 
Estimates 


Gravity 

Field  Estimation 


Convergence? . 


Adjustments  to 
Parameters 


Estimated  Htj 
Model  J 


Covariance 

Matrix 


Analysis  w/ 
PEP  &  Auxiliary 
Software 


Model  Analysis 


Figure  6.1-1:  Flow  Diagram  for  Gravity  Sensing  Mission  Simulations 

and  partial  derivatives,  the  parameter  estimation  module  formed  and  solved 
the  normal  equations  to  calculate  the  adjustments  to  the  parameters,  A  a. 
The  size  of  the  parameter  adjustments  and  their  uncertainties  determined  if 
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the  process  had  converged.  If  the  process  had  not  converged,  the  new 
parameters  were  used  to  propagate  the  next  iteration's  orbit  and  the  process 
was  repeated.  The  maximum  likelihood  parameter  estimates  and  the 
covariance  matrix  (A*1)  were  written  to  output  files  when  the  process  finally 
converged.  These  files  were  then  used  to  analyze  the  estimated  gravitational 
model  and  its  ability  to  predict  future  lunar  missions. 

Although  this  process  is  straight-forward,  it  is  a  very  time  and 
computer  memory  intensive  process.  A  typical  14  day  orbit  propagation  for 
an  8  X  8  spherical  harmonic  gravity  field  with  partials  calculated  for  79 
parameters  typically  took  over  eight  hours  on  a  Sim  Sparcstation  IPC.  Also, 
due  to  the  high  correlations  between  gravitational  parameters,  it  was  not 
unusual  for  a  run  with  different  "truth"  and  fit  models  to  require  at  least 
fifteen  iterations  to  converge  upon  a  solution.  The  convergence  criterion, 
perhaps  unnecessarily  stringent,  required  the  adjustments  to  be  less  than  1% 
of  the  parameter  uncertainties. 


6.2  Fifth  Degree  Harmonic  Truth  and  Fit  Model  Test  Case 

Since  PEP  was  modified  to  include  lunar  mascons  and  expand  the 
degree  and  order  spherical  harmonic  expansions  allowable,  several  test  cases 
were  rim  to  ensure  that  PEP-D's  modules  were  operating  correctly.  The  first 
set  of  tests  involved  the  attempt  to  estimate  a  5  X  5  spherical  harmonic 
expansion  fit  model  for  observations  generated  with  a  5  X  5  spherical 
harmonic  expansion  truth  model.  Since  the  degree  and  order  of  the 
estimated  fit  and  truth  models  were  the  same,  this  test  was  expected  to  be  a 
good  warm-up  exercise  for  future  runs. 

For  these  tests,  a  single  lunar  polar  200  km  altitude  near-circular 
satellite  orbit  was  numerically  integrated  for  twenty-eight  days.  The  initial 
selenographic  orbital  elements  are  listed  in  Table  6.2-1,  where  elements  in 
parentheses  are  in  PEP's  internal  units  (angles  referred  to  the  mean  equinox 
and  equator  of  the  earth  of  1950.0).  A  zero  degree  initial  mean  anomaly  was 
selected  first,  placing  the  satellite  directly  over  the  lunar  north  pole.  Due  to  a 
singularity  in  PEP's  algorithm  for  the  central  body  harmonic  effects,  it  could 
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not  propagate  this  orbit.  Because  of  perturbation  forces,  it  would  be  highly 
unlikely  that  this  condition  would  be  repeated  during  the  middle  of  an  orbit 
propagation  The  mean  anomaly  was  adjusted  one  degree  and  the  process 
restarted  with  no  further  difficulties. 

Table  6.2-1:  Satellite  Initial  Conditions  for  5  X  5  Spherical  Harmonic  Test  Case 


ao 

1938  km  (1.295472995  x  IO*5  AU) 

eo 

0.05 

io 

90°  (103.1048493849350°) 

Qo 

90°  (304.1996805394721°) 

ao 

90°  (72.2394269798987°) 

Mo 

1° 

The  starting  time  for  the  orbit  integration  was  16  May  1968  or  Julian 
Date  2,440,001.5  0  hour  Coordinate  Time.  This  initial  epoch  was  selected  over 
an  epoch  in  1996  or  1997,  when  a  lunar  gravitational  sensing  mission  might 
occur,  because  it  is  near  the  beginning  of  the  n-body  file  supplied  by  the  SAO. 
This  file  is  read  and  interpolated  from  during  the  integration  and  observation 
modules,  so  computer  time  was  saved  by  selecting  an  initial  epoch  near  the 
beginning  of  the  file.  Other  than  saving  computer  time,  the  choice  of  initial 
epoch  should  not  effect  file  simulations  of  this  thesis. 

The  lunar  satellite's  orbit  was  propagated  using  the  first  five  degree  and 
order  coefficients  from  the  1980  Bills  and  Ferrari  16  X  16  lunar  spherical 
harmonic  model  [12].  These  coefficients  are  given  in  Table  6.2-2.  The  central 
body  attraction  of  the  moon  and  the  perturbing  attractions  of  the  earth  and 
sun  were  included  in  the  Adams-Moulton  numerical  integration  with  a  step 
size  of  2*14  day.  This  allowed  approximately  1,451  steps  per  2.125504  hour 
orbit.  This  was  an  extremely  small  step  size  to  characterize  fifth  degree  and 
order  harmonics  and  ensure  numerical  stability  in  the  integration  of  motion. 

The  observation  module  then  created  Doppler  observations  of  the 
satellite  from  NASA  DSN's  Goldstone  site  every  minute.  Four  different 
observation  scenarios  were  generated  for  this  test  case.  Two  involved  range 
and  Doppler  observations  and  two  involved  Doppler  observations  only.  In 
two  cases  observations  were  occulted  by  the  moon  and  in  the  other  two  cases 
the  observations  continued  through  a  fictitious  "transparent"  moon. 


96 


Table  62-2:  5X5  Spherical  Harmonic  Truth  Model  Coefficients 


Harmonic 

x  10-6 

Harmonic 

xlO*6 

J2 

202.431 

■I 

C21 

-0.07 

-0.00 

C22 

34.49 

mm 

0.03 

h 

8.8897 

C31 

21.96 

S31 

6.63 

C32 

14.14 

S32 

4.76 

C33 

15.87 

S33 

-2.45 

J4 

-11.73 

C41 

-4.82 

S41 

C42 

-8.13 

S42 

Uil]  11  I11  111 1  11 

C43 

0.48 

S43 

-14.43 

C44 

-3.50 

S44 

-0.55 

Js 

2.388 

C51 

-9.66 

S51 

-1.53 

C52 

3.71 

S52 

-2.35 

C53 

-0.39 

S53 

4.91 

C54 

0.56 

S54 

-6.58 

C55 

-6.69 

S55 

11.60 

The  initial  test  cases  recreated  the  measurements  of  the  Lunar  Orbiter 
missions  (Section  1.3).  The  orbits  were  selected  to  simulate  orbits  proposed 
for  the  Lunar  Observer  or  Lunar  Scout  missions  [16,  25,  39].  The  observation 
scheme,  duplicating  the  processing  of  Lunar  Orbiter  data,  used  the  DSN  S- 
band  frequency  (2.115  GHz)  for  collecting  Doppler  range  rate  observations  for 
sixty-second  intervals.  No  transponder  frequency  translation  was  simulated. 
These  observations  had  a  quoted  accuracy  of  1  mm/sec  [29],  which  is  the 
accuracy  achievable  for  integer  Doppler  counts  over  the  sixty-second  interval. 
Range  observations  were  also  generated  in  the  "truth"  model  observation 
file.  When  fitting  to  only  Doppler  observations,  the  range  observations  were 
given  negative  error  weights  so  that  they  would  be  ignored  in  the  estimation 
process.  The  two  way  range  measurement  was  included  on  the  same  signal 
used  for  the  Doppler  observable  and  had  a  three  meter  accuracy.  This 
accuracy  is  based  upon  the  ability  to  detect  a  20  nanosecond  time  delay  for  a 
two-way  phase  shift  keyed  code. 


97 


LUNAR  GRAVITATIONAL  FIELD  ESTIMATION  AND  SATELLITE  ORBIT  PREDICTION 


Over  the  twenty-eight  day  observation  period,  the  satellite  would  orbit 
over  the  entire  lunar  surface  as  the  moon  rotated  under  its  orbital  plane.  Due 
to  the  extremely  conservative  integration  step  size,  the  twenty-eight  day 
viewing  period,  and  the  calculation  of  partials,  these  initial  test  cases  quickly 
filled  up  the  available  disk  space,  resulting  in  a  system  crash.  Reevaluating 
the  simulation  requirements,  the  observation  period  was  reduced  to  fourteen 
days,  which  still  provides  complete  lunar  coverage  in  either  an  ascending  or 
descending  pass. 

After  the  "true"  observations  were  generated,  the  satellite  initial 
conditions  and  thirty-two  lunar  harmonic  coefficients  were  perturbed  to 
provide  the  first  guess  for  the  estimation  process.  Each  gravitational 
harmonic  coefficient's  absolute  value  was  increased  by  1.0  x  10*7.  The  orbital 
initial  conditions  were  altered  by  the  values  given  in  Table  6.2-3.  The  satellite 
orbit  was  then  numerically  integrated  along  with  partial  derivatives. 

Table  6.2-3:  Perturbed  Satellite  Initial  Conditions  for  Estimation  Model 


8a 

149.6  km  (1.0x10*”  AU) 

8e 

1.0  xlO*5 

Si 

-1.0  x  1(H  0 

8Q 

-1.0x10*4° 

S(o 

-1.0  x  10*4  ° 

8Mq 

1.0  x  10*4  ° 

These  tests  were  run  using  saved  partial  derivatives  of  motion.  Using 
this  feature,  partial  derivatives  were  calculated  during  the  first  iteration  step 
and  these  saved  partials  of  motion  were  used  to  determine  the  partial 
derivatives  of  observations  on  subsequent  iterations.  Computation  time  was 
saved  by  not  recalculating  the  partials  on  each  iteration.  This  approach 
worked  well  for  test  cases  with  the  same  "truth"  and  fit  models,  but  was  not 
effective  on  other  simulations.  In  general,  this  feature  is  very  effective  when 
the  process  is  close  to  convergence. 

The  fitting  process  revealed  that  all  of  the  harmonic  terms  included  in 
the  equations  of  motion  have  to  be  included  in  the  partial  derivative 
calculations.  By  default,  PEP  was  only  calculating  the  second  degree  harmonic 
effects  upon  the  partials.  This  provides  sufficient  accuracy  for  earth  satellites 
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where  the  J2  effect  dominates  motion/  but  it  was  not  sufficient  for  lunar 
satellites.  Additionally/  the  first  scenario's  slow  convergence  raised  concerns 
about  the  choice  of  initial  osculating  orbital  elements.  For  a  circular  orbit,  the 
argument  of  perilune  (co)  is  undefined.  For  the  near-circular  orbit,  the 
perilune  position  was  very  difficult  to  determine.  The  argument  of  perilune 
and  mean  anomaly  were  difficult  to  estimate  because  of  their  high 
correlation.  PEP  software  does  not  currently  use  equinoctal  orbital  elements 
[e  sin(Q  +  to),  e  cos(£2  +  to)],  but  it  provides  the  option  of  replacing  the  Q,  to,  and 
M  orbital  elements  with  the  sum  of  the  angle  elements  Q,  Q  +  co,  £2  +  co  +  M. 
After  the  initial  orbital  elements  were  converted  to  this  form,  their 
correlations  were  reduced  and  the  parameter  estimation  routine  had  no 
difficulty  estimating  these  parameters. 

After  overcoming  the  previously  mentioned  difficulties,  the 
simulations  converged  in  three  iterations  for  all  four  scenarios.  After  these 
test  cases,  a  fifth  scenario  was  nm  in  which  radar  biases  for  the  time  delay  and 
Doppler  shift  observations  were  estimated  along  with  the  gravitational 
parameters  and  satellite  initial  conditions.  This  scenario  evaluated  the 
estimation  procedure's  ability  to  handle  measurement  bias  estimates  and  also 
converged  in  three  iterations. 


6.2.1  Scenario  One:  No  Occupation,  Range  and  Range  Rate 

The  first  scenario  calculated  a  spherical  harmonic  fit  model  based  on 
20,156  time  delay  (range)  and  20,156  Doppler  shift  (range  rate)  observations 
over  the  fourteen  day  period.  Observations  were  processed  every  sixty 
seconds  for  the  entire  period  since  there  were  no  horizon  or  occultation 
restrictions.  The  thirty-two  gravitational  harmonic  coefficients  and  six  initial 
conditions  were  all  estimated  within  approximately  ten  digits  of  accuracy  for 
each  parameter  (six  digits  for  angular  initial  conditions).  The  observation 
residuals  divided  by  the  assumed  measurement  errors,  referred  to  as  the  non- 
dimensional  "fit  residuals",  had  a  root  mean  square  (rms)  value  of  7.45263  x 
1(H,  essentially  zero.  The  theoretical  rms  for  these  cases  is  zero,  as  opposed  to 
the  value  of  one  mentioned  in  Section  5.2  because  the  "true"  observations 
were  exact  (without  measurement  noise)  and  the  truth  and  fit  models  are  of 
the  same  degree  and  order.  The  nm  also  calculated  correlations  as  high  as 
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-0.998  between  the  C33  and  C53  coefficients  and  $32  and  S52  coefficients. 
Fifteen  of  the  703  correlations  were  greater  in  magnitude  than  0.95.  Initially, 
these  high  parameter  correlations  seemed  as  if  they  were  due  to  a  resonance 
in  a  particular  harmonic  frequency,  but  later  fits  to  higher  degree  models  also 
demonstrated  high  correlations  for  the  higher  degree  terms. 

6.2.2  Scenario  Two:  No  Occupation,  Range  Rate 

This  scenario  excluded  the  20,156  time  delay  (range)  observations  from 
the  previous  rim's  fourteen  day  period  in  its  estimation  of  the  spherical 
harmonic  fit  model.  Once  again  the  thirty-eight  parameters  were  all 
estimated  with  approximately  ten  digits  of  accuracy  and  the  fit  residuals  were 
essentially  zero,  with  an  rms  value  of  8.83418  x  10*4.  The  highest  parameter 
correlation  was  again  -0.998  between  C33  and  C53.  The  S33S53  and  S32S52 
correlations  were  the  second  highest  at  -0.997.  The  S33  and  S53  correlation 
was  up  from  -0.996  in  the  previous  scenario.  Once  again  fifteen  of  the  703 
correlations  were  greater  in  magnitude  than  0.95. 

6.23  Scenario  Three:  Occupation,  Range  and  Range  Rate 

After  the  previous  scenarios  were  completed,  the  "truth"  model 
observations  were  recreated  with  lunar  occupations.  This  third  scenario  then 
calculated  its  fit  model  based  upon  15,672  time  delay  (range)  and  15,672 
Doppler  shift  (range  rate)  observations  over  the  fourteen  day  period.  Lunar 
occupations  eliminated  3  days,  2  hours  and  44  minutes  of  observations  over 
the  fourteen  day  period.  Again  the  thirty-eight  parameters  were  estimated  to 
within  ten  digits  accuracy  of  their  "true"  values.  With  the  reduced 
observability,  fit  residuals  increased,  having  an  rms  value  of  2.54901  x  10'3. 
Although  these  residuals  are  larger,  the  fit  model  is  still  essentially  exact.  The 
parameter  correlations  for  this  scenario  essentially  repeat  those  of  the  first 
scenario  with  the  C33C53  and  5 32 §52  coefficients  having  the  highest 
correlation  (-0.998).  As  in  both  previous  cases,  fifteen  of  the  703  correlations 
were  greater  than  0.95. 
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6.2.4  Scenario  Four:  Occulta  tion,  Range  Rate 

This  scenario  excluded  the  15/672  time  delay  (range)  observations  from 
the  previous  run  in  its  estimation  of  the  gravitational  coefficients  and  orbital 
initial  conditions.  The  scenario's  fit  model  agreed  with  the  truth  model  to 
about  nine  places.  The  fit  residuals  for  this  scenario  were  up  from  the 
previous  scenario  with  an  rms  value  of  6.86066  x  10*3.  The  highest  parameter 
correlation  was  the  C33C53  correlation  (-0.999).  The  532352  and  333353 
correlations  were  once  again  close  behind  (-0.998).  As  before,  fifteen  of  the  703 
correlations  were  greater  than  0.95. 


6.2.5  Scenario  Five:  Occultation,  Range  and  Range  Rate  with  Biases 

This  scenario  used  the  15,672  time  delay  (range)  and  15,672  Doppler 
shift  (range  rate)  observations  from  the  third  scenario  and  attempted  to 
estimate  the  thirty-two  gravitational  harmonic  coefficients  and  six  orbital 
initial  conditions  as  well  as  measurement  biases  for  the  time  delay  and 
Doppler  shift  observations.  This  scenario  estimated  the  gravitational 
coefficients  and  initial  conditions  with  about  nine  places  of  accuracy. 
Additionally  it  estimated  a  1.43121  x  10*9  second  time  delay  bias  and  a  -4.37639 
x  10"9  Hz  Doppler  shift  bias.  No  biases  were  included  in  the  observations  an 
this  run  essentially  estimated  zero  biases.  The  fit  residuals  were  lower  than 
the  fourth  scenario's  with  an  rms  value  of  6.10319  x  10-3,  also  essentially  zero. 
The  C33  and  €53  coefficients  once  again  had  a  -0.999  correlation  with  the 
§32^52  and  333353  parameters  yielding  the  second  highest  correlations  (-0.998). 
Fifteen  of  the  780  correlations  were  greater  than  0.95. 


63  Earth-Based  Doppler  Observable  Mascon  Test  Cases 

The  second  set  of  tests  included  two  small  mascons  in  the  truth  model 
and  increased  the  degree  of  the  spherical  harmonic  fit  model.  This  test  case 
was  run  to  determine  whether  PEP-D  could  estimate  a  spherical  harmonic 
model  based  on  observations  with  a  different  truth  model  and  to  determine 
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whether  the  sensing  geometry  affected  the  estimate's  ability  to  model  the 
lunar  gravitational  held. 

For  the  first  scenario,  the  previous  5x5  spherical  harmonic  expansion 
(Table  6.2-2 )  was  augmented  with  two  mascons  of  equal  strength  placed  at  the 
limbs  of  the  moon  on  the  equator.  For  the  second  scenario,  the  5  X  5  model 
was  augmented  with  mascons  of  equal  strength  placed  on  the  front  and  back 
faces  of  the  moon,  still  at  the  equator.  Mascons  were  included  in  pairs  to 
preserve  the  center  of  mass  of  the  moon  and  avoid  solving  for  first  degree 
harmonic  coefficients.  The  strongest  surface  disk  estimated  in  the  Wong 
model  [47]  was  selected  for  the  strength  of  these  mascons  whose  strengths  and 
locations  are  listed  in  Table  6.3-1. 


Table  63-1:  Mascon  Test  Case  -  Masccm  Placement 


Mascon  on  Limb 

Mascon  on  Face  j 

1 

2 

1 

2 

ESI 

7.212  x  10-6 

7.212  x  10-6 

7.212  x  10-6 

7.212  x  10-6 

_ Ei _ 

1680  km 

1680  km 

1680  km 

1680  km 

!E9 

90.0° 

180.0° 

I 

The  same  lunar  polar  200  km  altitude  near-circular  satellite  orbit  was 
numerically  integrated  for  fourteen  days  over  these  two  lunar  gravitational 
fields.  For  both  cases,  the  orbit  integration  was  once  again  started  on  16  May 
1968  and  used  the  same  satellite  initial  conditions  as  in  the  previous  tests 
(Table  6.2-1).  The  lunar  polar  satellite's  motion  was  detected  by  the  same 
earth-based  Doppler  shifted  observables  used  for  the  previous  tests. 
Observations  were  again  interrupted  by  lunar  occultations.  To  reduce  the  size 
of  the  computer  files  generated  in  the  simulation,  the  Doppler  count  interval 
was  changed  from  sixty  seconds  to  one  hundred  and  twenty  seconds. 

These  tests  were  run  to  determine  whether  earth  observations  are 
sufficient  to  estimate  the  lunar  gravitational  field  despite  the  mascons' 
location.  Mascons  cause  radial  disturbing  accelerations  as  satellites  pass  over 
them  and  Doppler  observations  sense  relative  velocity  along  the  line-of-sight. 
Doppler  observations  should  therefore  have  an  easier  time  sensing  mascon 
disturbances  when  the  observing  site,  observed  body,  and  mascon  are  all 
aligned.  These  two  lunar  mascon  scenarios  present  two  geometrical  extremes 
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for  Doppler  sensing.  In  the  first  case,  the  disturbing  acceleration  due  to 
passage  over  a  masoon  on  the  limb  is  orthogonal  to  the  Doppler  line-of-sight. 
In  the  second  case,  the  disturbing  acceleration  and  Doppler  line-of-sight  are 
aligned  as  the  orbiting  satellite  passes  over  the  near  face's  mascon. 

These  test  cases  were  initially  run  using  the  saved  partials  feature  of 
PGP  mentioned  in  Section  6.2.  For  the  first  scenario  full  partials  were 
recalculated  on  the  fifth  iteration  and  the  parameters  converged  on  the  ninth 
iteration.  The  parameters  from  this  solution  were  then  run  with  full  partials 
and  they  required  six  more  iterations  to  converge.  For  the  second  scenario, 
the  first  five  iterations  were  run  using  the  saved  partials  feature  since  the 
"semi-convergence"  of  the  previous  case  had  not  been  discovered.  After  the 
fifth  iteration,  the  saved  partials  feature  was  abandoned  and  this  scenario  also 
required  fifteen  total  iterations  to  converge. 

The  mascon  test  cases,  requiring  fifteen  iterations  to  converge, 
demonstrated  the  difficulty  of  estimating  gravitational  parameters  when  the 
"truth"  and  fit  models  differed.  Additionally,  these  scenarios  revealed  the 
limitations  of  PEP'S  saved  partials  feature.  Since  it  was  not  clear  that  saved 
partials  were  helping  the  estimation  process,  the  method  was  abandoned. 
This  significantly  increased  the  amount  of  computer  time  and  the  size  of  the 
memory  files  required  for  the  simulation  process. 


6.3.1  Scenario  One:  Limb  Mascons 

As  mentioned  previously,  this  scenario  converged  in  fifteen  iterations 
to  a  degree  and  order  eight  spherical  harmonic  model.  Appendix  E's  Table  E-l 
lists  the  estimated  harmonic  coefficients  for  this  fit  model.  Additionally  this 
run  estimated  the  satellite  osculating  orbital  initial  conditions.  Table  6.3.1-1 
compares  the  true  initial  conditions  with  their  estimated  values. 

The  fit  residuals  (observation  residuals  divided  by  assumed 
measurement  noise)  for  this  model  have  an  rms  of  7.36725.  For  the  first  time 
this  non-dimensional  statistic,  which  measures  how  well  the  estimated 
model  agrees  with  the  observed  behavior,  is  greater  than  one.  This  shows 
that  this  spherical  harmonic  model  is  inefficient  at  modeling  a  gravitational 
field  with  mascon  anomalies.  An  analysis  of  the  parameter  correlations 
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reveals  that  four  pairs  of  correlations  have  cross  correlations  greater  than 
-0.999:  §66586/  ^66^86/  ^55^75/  and  §55§75-  The  <£54^74  and  §54§74 

correlations  also  had  an  extremely  high  correlation  of  -0.989.  Twenty-two  of 
the  3,403  correlations  were  greater  in  magnitude  than  0.90  with  nine  of  these 
greater  than  0.95.  These  correlations  for  the  same  200  km  lunar  polar  orbiter 
provided  an  indication  that  the  high  parameter  correlations  discovered  in 
Section  6.2  were  not  due  to  the  resonance  of  a  particular  harmonic  frequency. 


Table  6.3.1-1:  Limb  Mascon  Orbital  Initial  Condition  Estimates 


"True" 

Fit 

A 

ao 

1938.0  km 

1938.1106  km 

110.6  meters 

eo 

0.05 

0.050091 

0.000091 

io 

103.1048° 

103.1067° 

0.0019° 

Oo 

304.1997° 

304.1939° 

-0.0058° 

(Qhd) o 

16.4391° 

16.6468° 

0.2077° 

(D-mih-M)o 

17.4391° 

17.6389° 

0.1998° 

6.3.2  Scenario  Two:  Face  Mascons 

This  scenario  fit  the  observed  behavior  to  an  8  X  8  spherical  harmonic 
expansion  in  fifteen  iterations.  Table  E-2  in  Appendix  E  lists  the  estimated 
harmonic  coefficients  for  this  fit  model.  This  run  also  estimated  the  satellite 
orbital  initial  conditions  for  the  observed  lunar  satellite.  .Table  6.3.2-1 
compares  the  true  initial  conditions  with  the  estimated  initial  conditions. 

This  case's  non-dimensional  fit  residuals  had  a  root  mean  square  of 
11.8025,  once  again  an  order  of  magnitude  greater  than  desired  for  an 
estimated  model.  The  increase  in  the  residuals'  rms  over  the  previous 
scenario  is  most  likely  due  to  the  lack  of  observations  for  the  far-side  mascon. 
The  Doppler  observations  in  this  scenario  sense  the  radial  accelerations  due 
to  the  mascon  on  the  near  face  but  cannot  account  for  the  lunar  far-side 
perturbations  with  the  8X8  spherical  harmonic  expansion.  In  the  previous 
case,  since  the  mascon  disturbing  accelerations  were  orthogonal  to  the 
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Doppler  line-of-sight,  the  observations  did  not  sense  large  disturbances  and 
had  an  easier  time  fitting  the  small  disturbances  to  the  fit  model. 

An  analysis  of  the  parameter  correlations  once  again  reveals  that  four 
pairs  of  correlations  are  greater  than  -0.998:  366586/  ^66^86/  C  55^75/  and 
555375.  The  C54C74  and  354574  correlations  are  again  very  highly  correlated 
at  -0.989  and  -0.988  respectively.  Of  the  3,403  correlations,  nineteen  were 
greater  in  magnitude  than  0.90  and  nine  of  these  were  greater  than  0.95. 


Table  6.3.2-1:  Face  Mascon  Orbital  Initial  Condition  Estimates 


'True" 

Fit 

A 

ao 

1938.0  km 

1937.9781  km 

-21.9  meters 

eo 

0.05 

0.049956 

-0.000044 

io 

103.1047° 

103.1015° 

-0.0032° 

Do 

304.1997° 

304.2022° 

0.0025° 

(D+co)o 

16.4391° 

16.3311° 

-0.1080° 

17.4391° 

17.3624° 

-0.0767° 

6.4  Truth  Model  Development 

Since  the  true  lunar  gravitational  field  is  not  precisely  known,  a  lunar 
gravitational  "truth"  model  was  developed  for  this  thesis.  This  truth  model 
was  then  used  to  evaluate  various  proposed  lunar  gravitational  sensing 
schemes.  Previous  lunar  missions  have  discovered  mascons  on  the  near  side 
of  the  moon  [35].  There  is  no  conclusive  evidence  regarding  mascons  on  the 
far  side  of  the  moon,  since  there  are  no  observations  of  a  satellite's  motion 
over  the  far  side.  For  this  thesis'  truth  model,  mascons  were  placed  on  both 
the  near  and  far  sides.  These  mascons  augmented  the  previously  used  5X5 
spherical  harmonic  expansion  (Table  6.2-2). 

From  their  surface  layer  representation  of  the  lunar  gravitational  field, 
Wong  et.  al.  identified  major  lunar  mass  anomalies  with  their  selenographic 
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features  [47  Table  2].  Wong's  group  modeled  the  lunar  gravitational  field 
with  50  km  radii  surface  disks  distributed  around  the  lunar  surface  (Section 
2.3.2).  Each  surface  feature  was  therefore  represented  by  several  disks 
depending  upon  the  size  of  the  feature.  The  seven  most  influential  mass 
anomalies,  requiring  forty-nine  surface  disks,  were  selected  for  the  truth 
model's  near-side  mascons.  Since  PEP-D  does  not  model  mascons  as  surface 
disks,  these  forty-nine  disks  were  converted  to  point  masses  placed  58  km 
below  the  lunar  surface.  The  stenographic  positions  and  strengths  of  the  600 
surface  disks  [47  Figure  4  and  Table  6]  were  then  correlated  with  topographic 
and  gravitational  maps  depicting  the  major  surface  features  [26]  to  determine 
the  forty-nine  point  masses  required  for  the  "truth"  model. 

Since  all  of  the  previously  identified  mascons  are  on  the  near  side, 
there  was  no  scientific  basis  for  establishing  lunar  far-side  mascons.  For  this 
thesis,  far-side  mascons  were  developed  to  meet  two  requirements.  First,  they 
should  be  difficult  to  detect  from  earth-based  observations,  and  secondly  they 
should  not  alter  the  lunar  center  of  mass.  Since  a  mascon  placed  on  the  back 
face  of  the  moon  would  be  the  most  difficult  to  sense,  the  first  far-side  mascon 
was  centered  at  180°  longitude  and  0°  latitude.  Masses  comparable  to  the 
near-side  mass  anomalies  were  selected  for  these  point  masses.  Since  lines  of 
longitude  are  more  widely  spaced  at  the  equator  and  nineteen  point  masses 
were  used,  this  mascon  is  the  largest  and  strongest  of  any  mass  anomaly  in 
the  truth  model. 

The  translation  of  the  lunar  center  of  mass  determined  the  placement 
of  the  second  far-side  mascon.  This  balancing  mascon  was  composed  of  seven 
point  masses  centered  on  214.5°  longitude  and  -48.0°  latitude.  After  these  two 
far-side  mascons  were  added  to  the  near-side  mascons,  the  center  of  mass  was 
still  askew,  so  an  additional  point  mass  added  to  the  model.  The  entire 
gravitational  truth  model  is  contained  in  Appendix  D. 

The  lunar  polar  200  km  altitude  orbit  used  in  the  two  previous  test 
cases  was  then  numerically  integrated  for  twenty-eight  days  with  this  lunar 
gravitational  truth  model.  The  satellite  orbit  was  then  converted  to 
selenographic  orbital  elements  and  plotted  to  determine  whether  the 
mascons  had  drastically  altered  the  lunar  gravitational  field.  The  plots 
revealed  that  the  mascons  affected  the  satellite  when  it  was  in  fhe  prime 
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meridian  plane/  seven  and  twenty-one  days  into  the  orbit  propagation.  At  this 
point/  the  instantaneous  semi-major  axis  increases  noticeably.  The  lunar 
altitude/  however,  never  drops  below  100  km  over  the  integration  period. 
The  lunar  altitude  for  this  orbiter  oscillates  between  100  -  300  km.  When 
numerically  integrated  without  mascons,  the  satellite  orbit  exhibits  the  same 
altitude  oscillations,  but  the  semi-major  axis  does  not  show  the  peaks 
observed  in  the  mascon  case.  The  inclusion  of  mascons  does  not  affect  the 
orbital  inclination  and  has  a  limited  effect  on  the  variation  of  eccentricity 
over  the  twenty-eight  day  period.  This  analysis  showed  that  the  addition  of 
mascons  to  the  lunar  gravitational  field  did  not  detrimentally  affect  the 
satellite  orbit. 


Table  6.4-1:  'Truth"  Model  Major  Mass  Anomalies 


Mascon  Strength  Longitude  Latitude  #of 


xlO** 

Lunar  Mass 


Point 

Masses 


328°  -350°  28°  -  48°  11 


go  .250  17o  .  340  12 


52°  -  63°  12°  -  23° 


27°  -  38°  -17°  -  -19° 


347°  -353°  7°  -  13° 


318°  -325°  -29° --22° 


82°  -  93°  -8°  -  3 


Sea  of  Rains 


Sea  of  Sereni 


Sea  of  Crises 


Sea  of  Nectar 


Sea  of  Moisture  6.0 


Subtotal  75.2 


Difficult  to  Obs  40.025  175°  -185°  -8°  -  8°  19 


Balancin 


Point  Mass 


27.157  212°  -217°  -51°  - -45° 


1.114 


217.6° 


-63.4 


Subtotal  68.296 


6.5  Single  Orbiter,  Earth-Based  Doppler  Sensing  Scheme 

The  first  sensing  scheme  evaluated  with  the  lunar  gravitational  truth 
model  was  the  earth-based  Doppler  observation  of  a  single  lunar  orbiter.  This 
scheme  was  selected  because  it  was  used  to  obtain  the  current  lunar  gravity 
field  models.  This  scheme's  ability  or  inability  to  estimate  the  "true"  lunar 
gravitational  field  may  provide  an  idea  of  how  well  current  models  represent 
the  real  lunar  gravity  field.  The  lunar  orbiter  was  placed  in  a  polar  near¬ 
circular  orbit  which  would  provide  total  lunar  surface  coverage  over  the 
fourteen  day  observation  period.  This  orbit  attempted  to  recreate  the  "gravity 
sensing"  satellite's  orbit  in  the  Lunar  Observer  or  Scout  missions  [16, 25, 39]. 

Since  the  acceleration  due  to  gravity  is  larger  for  low  altitude  orbits,  its 
disturbances  are  easier  to  sense  and  the  lowest  possible  orbital  altitude  was 
desired  for  this  mission.  Unfortunately  this  desire  conflicts  with  the  desire  to 
observe  undisturbed  motion  for  at  least  one  lunar  period.  Very  low  altitude 
orbits  typically  require  re-boosting  to  keep  them  at  a  safe  distance  above  the 
lunar  surface.  Re-boost  maneuvers,  however,  disrupt  the  estimation 
procedure  by  introducing  new  forces  on  the  orbiting  body  which  are  difficult 
to  include  in  the  estimation  process.  A  100  km  altitude  polar  orbiter  was 
numerically  integrated  for  twenty-eight  days.  Plots  of  its  selenographic  orbital 
elements  revealed  that  this  orbiter  barely  remained  above  the  lunar  surface 
for  the  twenty-eight  day  period.  Since  this  orbiter  would  require  re-boosting 
during  the  twenty-eight  days,  the  orbital  altitude  was  increased  by  100  km. 
The  200  km  altitude  polar  orbiter  maintained  an  altitude  of  100-300  km  over 
the  twenty-eight  day  period.  This  satellite  orbit  was  a  compromise  satisfying 
both  the  low-altitude  and  no  re-boost  requirements. 


6.5.1  Single  Orbiter  'Truth"  Model  Observations 

The  single  polar  200  km  altitude  lunar  orbiter  motion  was  numerically 
integrated  for  fourteen  days  from  16  May  1968  (JD  2,440,001.5)  with  the 
satellite  initial  conditions  given  below.  The  orbit  was  propagated  using  the 
Adams-Moulton  numerical  integration  technique  with  a  step  size  of  2*13  day. 
The  integration  step  size  was  relaxed  from  previous  cases  to  save  computer 
disk  space,  especially  when  partial  derivatives  are  calculated  at  each  step. 
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Relaxing  the  integration  step  size  still  provided  over  725  steps  per  2.125504 
hour  orbit,  more  than  required  for  numerical  stability  and  adequate  for  the 
characterization  of  an  eighth  degree  and  order  harmonic  fit  model. 

Over  this  fourteen  day  period,  Doppler  observations  of  the  lunar 
orbiter  were  processed  every  sixty  seconds,  except  during  lunar  occultations, 
from  the  DSN  Goldstone  station  using  the  S-Band  frequency  (2.115  GHz). 
Transponder  frequency  translation  was  not  simulated  and  no  horizon 
constraints  were  imposed.  A  14  mHz  accuracy  was  assumed  for  the  Doppler 
observations,  roughly  equivalent  to  the  quoted  1  mm/sec  range  rate  accuracy 
for  sixty  second  intervals  [29].  The  simulations  of  this  mission  produced 
15,695  Doppler  shift  observations  over  the  fourteen  day  period. 

Table  65.1-1:  Satellite  Initial  Conditions  for  "Truth"  Model  Observations 


ao 

1938  km  (1.295472995  x  10-5  AU) 

eo 

0.05 

io 

90°  (103.1048493849350°) 

Qo 

90°  (304.1996805394721°) 

■29 

90°  (72.2394269798987°) 

1° 

6.5.2  Eighth  Degree  and  Order  Fit 

The  first  estimation  run  attempted  to  fit  the  observations  to  an  8  x  8 
spherical  harmonic  expansion.  The  final  iteration  of  the  limb  mascon  test 
case  provided  the  initial  guess  of  parameters  and  initial  conditions  for  this 
case's  iterations.  This  saved  computer  time  since  both  cases  were  based  on 
observations  of  the  same  polar  orbiter  and  the  limb  mascon  case  had  already 
integrated  the  satellite's  motion  with  partial  derivatives. 

Due  to  the  significant  difference  between  the  model  used  to  generate 
the  "true"  observations  and  the  fit  model,  this  case  required  twenty-five 
iterations  to  converge  upon  a  solution.  The  estimated  harmonic  coefficients 
for  the  fit  model  are  listed  in  Table  E-3  in  Appendix  E.  The  simulation  also 
estimated  the  satellite  orbit's  initial  conditions  and  Table  6.5.2-1  compares  the 
true  initial  conditions  with  their  estimated  values. 
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The  observation  residuals  divided  by  the  assumed  measurement 
errors,  referred  to  as  the  non-dimensional  "fit  residuals",  had  a  root  mean 
square  value  of  365.954  for  this  model.  This  is  over  thirty  times  larger  than 
the  Face  Mascon  test  case  and  is  a  further  indication  of  the  8  X  8  spherical 
harmonic  expansion's  inability  to  fit  the  observed  satellite's  motion.  Eighteen 
of  the  3,403  parameter  correlations  were  greater  in  magnitude  than  0.90  and 
nine  of  these  were  greater  than  0.95.  Four  pairs  of  correlations  (C55C75, 
S 66^86/  ^ 66^86/  and  3 55 S 75)  were  greater  than  -0.998  and  two  pairs  (354S74  and 
C54C74)  were  greater  than  -0.988. 


Table  65.2-1:  8X8  Single  Orbiter  Initial  Condition  Estimates 


'True" 

Fit 

A 

ao 

1938.0  km 

1938.1154  km 

115.4  meters 

eo 

0.05 

0.047811 

-0.002189 

io 

103.1048° 

103.0998° 

-0.0050° 

Qo 

304.1997° 

304.1957° 

0.0040° 

(Q+cd)o 

16.4391° 

15.9600° 

-0.4791° 

(D-HiH-M)o 

17.4391° 

17.3912° 

-0.0479° 

6.5.3  Twelfth  Degree  and  Order  Fit 

The  estimated  gravity  model's  fit  residuals  show  that  the  8X8  spherical 
harmonic  expansion  was  not  a  close  fit  to  the  observed  behavior.  Mascons  in 
the  truth  model  result  in  very  high  frequency  gravitational  behavior  in  the 
local  area  (Section  2.3).  Since  the  spherical  harmonic  expansion  requires 
higher  degree  and  order  expansions  to  model  this  high  frequency  behavior,  a 
fit  to  a  12  X  12  spherical  harmonic  model  was  attempted. 

For  the  initial  parameter  guesses,  estimates  from  the  8  X  8  fit  model 
were  used  and  all  of  the  new  coefficients  (ninth  through  twelfth  degrees) 
were  set  to  zero.  The  satellite  initial  condition  estimates  for  the  8X8  fit,  used 
for  the  12  X  12  fit's  initial  iteration,  are  listed  in  Table  6.5.3-1.  The  same 
"truth"  model  observations  used  in  the  8  X  8  fit  were  used  for  this  fit. 
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Since  the  fit  model  was  increased  to  twelfth  degree  and  order,  the 
numerical  integration  step  size  was  reduced  to  2~14  days.  The  step  size  was 
changed  to  sample  the  higher  degree  harmonic  effects  more  often  with  1,451 
steps  per  2.125504  hour  orbit.  With  the  change  in  step  size  and  dramatically 
increased  number  of  estimation  parameters  and  parameter  partials,  the 
numerical  integration  required  approximately  23  hours  to  propagate  the 
satellite's  motion  and  partial  derivatives  for  fourteen  days. 

Table  65.3-1:  12  X  12  Fit  Initial  Guesses  for  Satellite  Initial  Conditions 


ao 

1.29555011876471  x  10-5  AU 

eo 

0.0478107828679364 

io 

103.099831198443° 

fio 

304.195650955435° 

(£2  +  cd)o 

15.9599887633931° 

(£2  +  co  +  M)o 

17.3911644150588° 

On  the  fitting  process'  second  iteration  several  of  the  parameter 
adjustments  were  an  order  of  magnitude  larger  than  their  previous  estimate. 
These  estimates  were  used  to  propagate  the  next  iteration's  orbit,  resulting  in 
increased  residuals  and  even  larger  parameter  adjustments.  By  the  fifth 
iteration,  the  normal  equations  could  not  be  solved  because  152  of  the  171 
diagonal  elements  in  the  coefficient  matrix  (A)  were  negative.  Since  the 
normal  equations  could  not  be  inverted,  the  fitting  process  was  abandoned. 


Table  65.3-2:  12  X  12  Spherical  Harmonic  Fit  Progression 


Iteration 

Pre-adjustment 
RMS  Residual 

Predicted  RMS 

Residual 

1 

365.954 

112.113 

2 

10,433.0 

99.5596 

3 

42,590.9 

192.212 

4 

106,490 

1,999.54 

5 

760,314,000 

16,484,700 

Table  6.5.3-2  shows  how  the  parameter  adjustments  kept  leading  the 
process  further  and  further  from  a  solution.  The  fitting  process  started  with 
observation  residuals  divided  by  the  assumed  measurement  errors  with  a 
root  mean  square  value  of  366.  After  solving  the  normal  equations,  PEP 
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predicted  a  non-dimensional  rms  residual  of  112.  The  predicted  residual  is 
the  computed  pre-adjustment  residual  less  the  sum  of  the  observable  partial 
derivatives  times  the  parameter  adjustments.  Instead  of  reducing  the 
residuals  by  a  third,  the  new  parameter  estimates  produced  residuals  with  a 
non-dimensional  rms  almost  thirty  times  greater  than  in  the  previous 
iteration!  Unfortunately  instead  of  making  it  easier  to  fit  to  the  observed 
lunar  orbiter  motion,  increasing  the  degree  and  order  of  the  spherical 
harmonic  model  made  it  more  difficult  to  solve  the  normal  equations  and 
converge  upon  a  fit  solution.  This  may  have  occurred  because  of  the  limited 
observability  of  lunar  far-side  motion  and  high  parameter  correlations. 


6.6  Tenth  Degree  Harmonic  Truth  and  Fit  Model  Test  Case 

After  the  previous  convergence  difficulties,  a  test  was  run  to  estimate  a 
10  X  10  spherical  harmonic  expansion  for  observations  generated  with  the 
same  degree  and  order  truth  model.  Since  the  previous  case  was  the  first 
attempt  to  fit  to  harmonic  expansions  higher  than  eighth  degree  and  order, 
this  test  would  verify  whether  or  not  the  difficulties  encountered  were  due  to 
the  estimation  program  or  the  spherical  harmonic  expansion's  ability  to  fit  to 
the  observed  behavior. 

The  same  lunar  polar  200  km  altitude  near-circular  orbit  used  in 
previous  runs  was  used  to  generate  Doppler  "truth"  observations  with  1 
mm/sec  accuracy.  This  orbit,  again  integrated  for  fourteen  days  from  16  May 
1968,  used  the  orbital  initial  conditions  listed  in  Table  6.2-1  and  a  10  X  10 
spherical  harmonic  lunar  gravity  truth  model.  This  truth  model  used  the 
1980  Bills  and  Ferrari  coefficients  [12]  for  the  initial  five  degrees  (Table  6-2.2). 
Gravitational  coefficients  developed  at  JPL  by  Alex  Konopliv  for  his  50  X  50 
spherical  harmonic  expansion  were  used  for  the  sixth  through  tenth  degree 
harmonic  coefficients  [27],  because  the  Bills  and  Ferrari  coefficients  [12]  were 
not  manually  entered  into  the  PEP  input  stream.  The  coefficients  used  are 
listed  in  Table  E-4  in  Appendix  E. 

After  generating  the  "truth"  model  observations,  the  satellite  initial 
conditions  and  one  hundred  and  seventeen  lunar  harmonic  coefficients  were 
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perturbed  to  provide  the  first  guess  for  the  estimation  process.  As  in  the 
previous  fit  to  a  pure  harmonic  truth  test  case,  each  gravitational  harmonic 
coefficient's  absolute  value  was  increased  by  1.0  x  10*7  and  the  orbital  initial 
conditions  were  perturbed  by  the  values  given  in  Table  6.2-3.  The  satellite 
orbit  was  then  numerically  integrated  with  partial  derivatives. 

In  the  first  iteration  the  pre-adjustment  observations  residuals  divided 
by  the  assumed  measurement  error  had  a  root  mean  square  value  of  715,121! 
As  a  result,  the  calculated  parameter  adjustments  were  orders  of  magnitude 
larger  than  the  guessed  parameters.  After  the  first  iteration,  the  lunar  orbiter 
eccentricity  grew  from  0.05001  to  0.24325  and  the  semi-major  axis  was 
increased  by  532  km.  These  initial  conditions  were  then  propagated  to  create 
the  next  iteration's  observations.  In  the  second  iteration,  the  normal 
equations  could  not  be  solved  because  51  of  the  123  diagonal  elements  in  the 
coefficient  matrix  (A)  were  negative.  Once  again  since  the  normal  equations 
could  not  be  inverted,  the  estimation  process  was  abandoned.  Rather  than 
converging  upon  a  solution.  Table  6.6-1  demonstrates  how  quickly  the  process 
diverged. 


Table  6.6-1:  10  X 10  Spherical  Harmonic  Fit  Progression 


Iteration 

Pre-adjustment 
RMS  Residual 

Predicted  RMS 

Residual 

1 

715,121 

248,467 

2 

3.20623  x  1010 

2.11726  xlO7 

Since  the  previous  test  case  (Section  6.2)  encountered  difficulty  when 
an  incomplete  set  of  harmonic  coefficients  were  used  to  generate  partial 
derivatives  of  the  satellite  motion,  this  cases'  observation  partial  derivatives 
were  checked  using  a  finite  difference  method  (Section  3.4.4).  This  check 
verified  that  the  observation  partials  were  correct. 

Since  the  estimation  software  was  operating  properly,  the  gravitational 
model's  sensitivity  to  initial  guesses  seemed  to  cause  tne  convergence 
difficulties.  Apparently  the  initial  guesses  resulted  in  theoretical  observation 
values  so  different  from  the  observed  behavior  that  correct  parameter 
adjustments  could  not  be  determined.  Since  the  estimation  of  higher  degree 
and  order  expansion  fit  models  is  extremely  sensitive  to  die  initial  parameter 
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guesses,  fit  models  above  eighth  degree  and  order  were  not  attempted  for  the 
remainder  of  this  thesis'  simulations. 

Smaller  parameter  adjustments  from  iteration  to  iteration  could  help 
the  process  converge.  Incorporating  a  priori  estimates  of  the  lunar 
gravitational  field  and  their  uncertainties  would  reduce  the  size  of  the 
parameter  adjustments  (see  Section  5.1.4).  Smaller  parameter  adjustments 
could  also  be  used  by  under-weighting  the  calculated  adjustments,  i.e.  only 
adjusting  the  initial  guesses  by  2/3  of  the  calculated  adjustment  on  the  first 
iteration,  etc.  Convergence  could  also  be  aided  by  first  estimating  an  8  X  8, 
then  a  9  X  9,  and  then  a  10  X  10  degree  and  order  fit  model  to  the  observed 
behavior,  using  the  estimates  from  one  model  as  initial  guesses  for  the  next 
higher  degree  model. 


6.7  Dual  Orbiter,  Bent  Pipe  Doppler  Sensing  Scheme 

The  next  sensing  scheme  evaluated  with  the  lunar  gravitational 
"truth"  model  featured  two  lunar  satellites  and  the  simulation  of  a  bent  pipe 
Doppler  sensing  scheme  between  earth-based  sites  and  the  two  satellites.  This 
sensing  scheme  was  selected  since  it  is  one  of  the  schemes  being  considered  by 
NASA  for  future  lunar  gravitational  sensing  missions. 

This  sensing  scheme  uses  a  low  altitude  circular  polar  "gravity 
sensing"  satellite.  A  coplanar  elliptical  "viewing"  satellite  makes  lunar  far- 
side  observations  of  the  "gravity  sensing"  satellite's  motion  [16,  39,  40].  The 
polar  200  km  altitude  near-circular  lunar  orbiter  used  in  previous 
simulations  is  once  again  the  "gravity  sensing"  satellite  for  this  dual  orbiter 
sensing  scheme.  The  satellite  initial  conditions  and  numerical  integration 
parameters  for  this  simulation  were  the  same  as  for  the  single  orbiter,  earth- 
based  Doppler  observation  sensing  scheme  and  are  given  in  Section  6.5.  The 
"viewing"  satellite  was  placed  in  a  450  km  X  7,000  km  altitude  elliptical  orbit 
with  a  10.06  hour  period.  This  satellite  was  given  the  same  initial  inclination 
(i)  and  longitude  of  the  ascending  node  (Q)  as  the  circular  satellite,  placing 
them  in  the  same  orbital  plane,  orbiting  the  moon  in  the  same  direction. 
With  this  selection  of  ascending  node,  the  far  side  of  the  moon  rotates 
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underneath  the  elliptical  orbit's  apolune  during  the  first  fourteen  days  of  orbit 
propagation.  This  maximizes  the  satellite's  ability  to  view  the  low  altitude 
arbiter's  far-side  motion. 

The  argument  of  peril une  was  selected  to  skew  the  ellipse  and  bring 
the  orbit's  apolune  position  out  of  the  earth  occupation  zone  (Figure  1.4-1). 
Keeping  the  satellite's  apolune  passage  out  of  the  occultation  zone  increases 
the  time  earth-based  sites  can  view  the  elliptical  satellite  for  bent-pipe 
measurements  as  well  as  data  transfers.  The  long  term  behavior  of  the 
osculating  orbital  elements  was  studied  to  determine  this  skew  direction. 
Formulas  for  the  doubly  averaged  effect  of  the  earth  upon  a  lunar  orbit 
revealed  that  if  the  initial  argument  of  peril  une,  co,  is  in  the  second  or  fourth 
quadrants,  the  polar  orbit's  eccentricity  will  decrease  and  <d  will  drift  to  the 
first  or  third  quadrant  [9].  Based  upon  this  analysis  and  die  desire  to  keep  the 
apolune  over  the  backside  of  the  moon,  a  perilune  angle  in  the  fourth 
quadrant  was  selected.  The  initial  selenographic  orbital  elements  for  this 
elliptical  satellite  are  given  in  Table  6.7-1,  with  the  angle  values  in 
parentheses  referred  to  the  mean  equinox  and  equator  of  the  earth  of  1950.0. 
These  initial  conditions  were  numerically  integrated  for  fourteen  days  from 
16  May  1968  (JD  2,440,001.5)  using  the  Nordsieck  variable  step  size  integration 
technique  (Section  3.1).  Since  the  elliptical  orbit  was  newly  propagated,  its 
numerical  integration  file  was  converted  to  selenographic  orbital  elements. 
Plots  of  these  elements  revealed  that  die  orbit  was  stable  and,  as  predicted, 
the  orbital  eccentricity  did  decrease  over  the  fourteen  day  period. 

Table  6.7-1:  Elliptical  Satellite  Initial  Conditions  for  "Truth"  Model  Observations 


ao 

5463  km  (3.651789462  x  IO"5  AU) 

eo 

0.05994874611 

io 

90°  (103.1048493849350°) 

Oo 

90°  (304.1996805394721°) 

COD 

315°  (297.2394269798987°) 

Mo 

1° 
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6.7.1  Dual  Orbiter  "Truth*  Model  Observations 

The  coherent  bent  pipe  observation  method  (Section  4.2.3)  was 
simulated  over  the  fourteen  days  the  satellites'  orbits  were  propagated.  One 
two-way  coherent  Doppler  loop  processed  observations  between  the  DSN 
Goldstone  site  and  the  elliptical  satellite.  Simultaneously,  the  elliptical 
satellite  generated  two-way  coherent  Doppler  observations  of  the  low-altitude 
orbiter.  Although  the  coherent  bent  pipe  scheme  would  be  interrupted  if  any 
of  the  links  were  occulted,  only  the  individual  links  in  this  simulation  were 
interrupted  by  lunar  occultations.  In  addition  to  the  bent  pipe  observation 
simulation,  the  DSN  Goldstone  site  generated  Doppler  observations  of  the 
near-circular  polar  satellite  during  its  near-side  passes.  All  of  these 
observations  were  simulated  using  the  2.115  GHz  S-band  frequency  with 
approximately  1  mm/sec  range  rate  accuracy  (14  mHz)  [29]. 

Simulating  this  sensing  scheme  produced  15,695  earth-based  Doppler 
observations  of  the  circular  satellite.  Because  the  elliptical  orbit  was  skewed, 
19,699  Doppler  observations  were  generated  between  it  and  the  DSN  station. 
Over  this  same  period,  10,680  Doppler  observations  were  simulated  between 
the  two  satellites.  Auxiliary  software  analyzed  these  observation  series  and 
their  occultation  periods  to  evaluate  the  far-side  lunar  coverage.  Over  the 
fourteen  days,  the  low  altitude  polar  orbiter  passed  behind  the  moon  114 
times.  For  most  of  these  occultations,  the  elliptical  satellite  viewed  this  far- 
side  passage.  The  moon  blocked  the  line-of-sight  between  the  two  satellites 
on  49  of  these  114  occasions  and  this  blockage  usually  only  affected  a  portion 
of  the  passage.  There  were  only  15  cases  in  which  the  line-of-sight  to  the 
elliptical  satellite  was  blocked  for  the  entire  far-side  passage.  Since  these  gaps 
in  observation  coverage  did  not  occur  for  sequential  far-side  passes,  the  46,074 
observations  should  provide  excellent  visibility  into  the  circular  orbiter's 
motion  on  the  lunar  far  side. 


6.7.2  Eighth  Degree  and  Order  Fit 

These  observations,  generated  with  the  lunar  gravitational  "truth" 
model,  were  then  fit  to  an  8  X  8  spherical  harmonic  expansion  representing 
the  moon's  gravitational  field.  The  initial  conditions  for  the  two  satellites 
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were  perturbed  from  their  true  values  for  the  initial  guess  in  the  fitting 
process.  Table  6.2-3  lists  the  perturbations  which  were  applied  to  both  the 
circular  and  elliptical  orbits'  initial  conditions.  The  true  values  of  the  5x5 
spherical  harmonic  expansion  used  in  the  truth  model  were  used  as  initial 
guesses  for  the  first  five  degrees  of  gravitational  coefficients  (Table  6.2-2). 
Gravitational  coefficients  from  Alex  Konopliv's  50  X  50  spherical  harmonic 
model  of  the  moon  [27]  were  used  for  the  sixth,  seventh,  and  eighth  degree 
coefficient  initial  guesses.  Table  E-4  in  Appendix  E  lists  these  coefficients. 

The  iteration  fit  did  not  converge  upon  a  solution  in  its  normal  sense. 
With  the  increased  observability  of  lunar  far-side  motion,  the  parameter 
estimation  routine  drove  the  non-dimensional  residuals  down  to  a  root 
mean  square  of  386.462  from  an  initial  23,763.6  in  the  first  four  iterations. 
Because  of  the  difference  between  the  "truth"  and  fit  models,  the  estimation 
routine  could  not  reduce  the  non-dimensional  residual  rms  below  360  as 
shown  in  Table  6.7.2-1.  Once  again,  the  predicted  residual  is  the  computed 
pre-adjustment  residual  minus  the  sum  of  the  observable  partial  derivatives 
times  the  parameter  adjustments. 


Table  6.72-1:  8x8  Spherical  Harmonic  Fit  Progression, 


Iteration 

Pre-adjustment 
RMS  Residual 

Predicted  RMS 

Residual 

1 

23,763.6 

817.104 

2 

2,983.65 

437.085 

3 

1,746.75 

383.599 

4 

386.462 

371.811 

5 

375.229 

366.58? 

6 

368.213 

363.927 

| 

| 

15 

361.210 

361.015 

j 

38 

361.099 

360.939 

Although  the  estimation  routine  had  minimized  the  fit  residuals,  its 
convergence  criteria  is  based  on  the  ratio  of  the  parameter  adjustments  to 
their  uncertainties  with  the  assumed  measurement  errors,  and  these  values 
remained  above  the  convergence  limit  of  0.01.  Since  the  convergence 
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criterion  did  not  signal  convergence,  the  estimation  iterations  continued  ad 
infinitum,  slowly  reducing  the  residual  rms  by  hundredths  and  thousandths. 
The  process  was  stopped  after  thirty-eight  iterations. 

The  estimated  harmonic  coefficients  for  the  fit  model  are  listed  in 
Table  E-5  in  Appendix  E.  The  simulation  also  estimated  the  orbital  initial 
conditions  for  both  of  the  satellites  in  the  sensing  scheme.  Tables  6.72-2  and 
6.72-3  compare  the  true  initial  conditions  with  their  estimated  values. 

The  estimated  model's  fit  residuals  had  a  root  mean  square  value  of 
361.099  as  listed  in  Table  6.7.2-1.  This  is  smaller  than  the  root  mean  square  of 
the  residuals  achieved  with  the  fit  to  a  single  lunar  polar  orbiter  (365.954). 
The  high  fit  residuals  are  an  indication  of  the  8  X  8  spherical  harmonic 
expansion's  inability  to  fit  the  observed  motion  due  to  the  "true" 
gravitational  field.  This  dual  satellite  observation  method  has  significantly 
reduced  some  of  the  high  parameter  correlations,  although  the  most 
significant  ones  still  remain.  For  this  fit  case  fifteen  of  the  3,916  parameter 
correlations  were  greater  in  magnitude  than  0.90  and  ten  of  these  were  greater 
than  0.95.  As  with  the  previous  8X8  spherical  harmonic  fit,  the  four  highest 
correlations  were  among  the  C55C75,  §55575,  §66§86/  and  ^66^86  pairs. 
Although  the  correlations  between  the  fifth  degree,  fifth  order  and  seventh 
degree,  fifth  order  terms  were  still  greater  than  -0.998,  the  other  correlations 
were  now  reduced  to  -0.995.  The  §54574  and  C54C74  pairs  were  again  the  next 
highest  correlated  (-0.988). 


Table  6.72-2:  Circular  Orbiter  Initial  Condition  Estimates 


'True" 

Fit 

A 

ao 

1938.0  km 

1937.9401  km 

-59.9  meters 

eo 

0.05 

0.049677 

-0.000323 

io 

103.1048° 

103.1033° 

-0.0015° 

Qo 

304.1997° 

304.1994° 

-0.0003° 

(Q+©)o 

16.4391° 

16.2649° 

-0.1742° 

(£2+o+M)o 

17.4391° 

17.4469° 

0.0078° 
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Table  6.72-3:  Elliptical  Orbiter  Initial  Condition  Estimates 


'True" 

Fit 

A 

ao 

5463.0  km 

5463.6291  km 

629.1  m 

eo 

0.599487 

0.599573 

0.000086 

k> 

103.1048° 

103.1068° 

0.0020° 

**> 

304.1997° 

304.1992° 

0.0005° 

COD 

2972394° 

2972542° 

0.0148° 

Mo 

O 

O 

r“4 

0.9997° 

-0.0003° 
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Chapter  Seven 

Lunar  Gravity  Field 
Estimation  Analyses 


7.1  Analysis  Techniques 

After  converging  upon  a  least  squares  maximum  likelihood  estimate 
of  gravitational  coefficients,  the  estimated  lunar  gravity  fields  were  analyzed 
to  compare  their  ability  to  estimate  the  "true"  gravity  field.  Since  spacebome 
navigation  depends  on  accurately  modeling  the  forces  acting  on  a  spacecraft, 
these  analyses  focus  on  the  effects  of  modeling  errors  on  lunar  navigation. 

The  root  mean  square  of  the  observation  residuals  does  not  provide  an 
adequate  analysis  of  the  global  lunar  gravity  errors  due  to  mismodeling,  since 
it  only  considers  the  areas  of  the  gravity  field  where  observations  were  made. 
This  underweights  the  far  side  of  the  moon  where  fewer  measurements  are 
taken.  In  Section  7.2  the  global  radial  accelerations  for  the  estimated  fit  lunar 
gravity  fields  are  compared  to  the  "true"  lunar  radial  accelerations.  The  limb 
mascon  and  face  mascon  estimated  fit  models  are  analyzed  and  the  two  8X8 
spherical  harmonic  expansion  fit  models  are  compared  to  the  lunar 
gravitational  "truth"  model  developed  in  Section  6.4. 

Next,  two  different  lunar  spacecraft  mission  phases  are  simulated  to 
evaluate  the  two  estimated  fits  to  the  lunar  gravitational  "truth"  model. 
These  analyses  show  how  state  errors  grow  using  the  estimated  gravity  field 
model  and  are  intended  to  simulate  the  real-world  consequences  of  planning 
and  executing  lunar  missions  with  a  mismodeled  gravity  field. 

In  Section  7.3  the  state  uncertainties  for  a  low  inclination,  low  altitude 
spacecraft  orbit  are  predicted  one  orbit  ahead  using  the  covariance  analysis 
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described  in  Section  5.2.  Stations  tracking  a  lunar  satellite  would  use  this 
analytical  method  during  far-side  passage  to  predict  the  spacecraft's  state 
uncertainties  upon  reemergence  from  the  backside.  Since  this  state 
uncertainty  prediction  is  based  upon  a  mismodeled  lunar  gravity  field,  the 
predicted  uncertainties  will  not  be  exact.  The  "true"  state  is  calculated  by 
numerically  integrating  the  spacecraft  orbit  using  the  lunar  gravity  "truth" 
model.  The  difference  between  the  spacecraft's  predicted  and  "true"  states  is 
the  error.  Comparing  the  "true"  state  with  the  predicted  states  reveals  the 
accuracy  of  the  estimated  lunar  gravity  field. 

A  lunar  landing  deorbit  maneuver  is  simulated  in  Section  7.4.  In  this 
analysis  the  deorbit  trajectory  from  a  low  inclination  circular  parking  orbit 
was  determined  from  the  estimated  lunar  gravity  field  models.  The  mission 
was  then  propagated  using  the  "true"  lunar  gravity  field.  Lunar  gravity 
mismodeling  resulted  in  a  spacecraft  position  error  by  the  time  the  spacecraft 
reached  the  lunar  altitude  for  Powered  Descent  Initiation  (PDI),  typically 
about  18  kilometers.  In  a  real  lunar  mission,  die  spacecraft  would  be  forced  to 
bum  extra  propellant  in  a  suboptimal  descent  trajectory  to  recover  from  these 
errors. 


Finally,  the  high  gravitational  coefficient  correlations  encountered  in 
the  estimation  process  are  analyzed  in  Section  7.5.  Three  different 
measurement  types  and  orbital  orientations  are  simulated  to  determine  their 
effect  on  the  correlations.  Reducing  the  coefficient  correlations  by 
introducing  new  observations  may  permit  the  process  to  converge  more 
quickly  and  estimate  a  more  accurate  gravitational  model.  Additionally,  the 
highest  correlations  are  studied  in  an  attempt  to  find  ways  to  allow  them  to  be 
estimated  independently. 


7.2  Global  Lunar  Radial  Acceleration  Analysis 

Although  the  estimation  process  provides  residual  statistics,  they  only 
provide  an  idea  of  how  well  the  estimate  fits  the  observations  since  the 
observations  are  not  available  over  the  entire  lunar  surface.  To  analyze  the 
fit  globally  a  software  program  was  written  to  calculate  the  radial  accelerations 
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for  a  gravity  field  model  everywhere  on  a  sphere  of  constant  radius.  Only  the 
radial  accelerations  due  to  gravity  harmonics  and/or  mascons  are  calculated 
in  this  program.  These  radial  accelerations  were  calculated  for  fit  and  truth 
models  for  all  of  the  cases  estimated  in  Chapter  Six  at  a  lunar  altitude  of  100 
km  for  grid  points  spaced  every  four  degrees  of  selenographic  latitude  and 
longitude.  Unfortunately,  since  the  commercial  graphics  package  was 
erroneously  plotting  lines  of  constant  radial  acceleration,  contour  plots 
comparing  the  fit  and  "truth"  models  are  not  currently  available. 

Based  on  the  calculated  radial  accelerations,  the  fit  and  "truth"  models 
were  compared  statistically.  The  radial  acceleration  errors  between  the 
"truth"  and  fit  models  were  calculated  for  each  grid  point.  The  root  mean 
square  error  between  the  two  models  was  then  calculated.  Because  of  the 
even  spacing  of  grid  points  in  selenographic  latitude  and  longitude  (chosen 
for  rectangular  contour  plots),  this  analysis  weights  radial  acceleration  errors 
more  heavily  in  the  polar  regions.  For  the  mascon  test  case  this  bias  does  not 
favor  either  scenario  since  the  disturbing  mascons  were  placed  on  the 
equator.  For  the  two  8X8  estimated  fit  models,  the  dual  orbiter  bent  pipe 
estimated  fit  might  have  an  advantage  since  the  "viewing"  satellite  observed 
the  "sensing"  satellite  as  it  passed  over  the  lunar  polar  regions.  Both  models, 
however,  were  trying  to  estimate  the  same  lunar  gravity  "truth"  model,  so  a 
comparison  of  the  rss  radial  acceleration  errors  is  still  valid. 


7.2.1  Limb  /  Face  Mascon  Analysis 

Based  on  the  global  lunar  radial  acceleration  errors,  the  face  mascon 
case  produced  a  slightly  better  estimated  gravity  field.  This  case  had  an  rms 
radial  acceleration  error  at  100  km  altitude  of  64.7804  milligals.  The  estimated 
fit  to  the  limb  mascon  case  had  radial  acceleration  errors  with  an  rms  of 
66.6636  milligals.  This  global  analysis  of  the  radial  accelerations  errors  reveals 
that  a  slightly  better  lunar  gravity  field  was  estimated  when  the  disturbing 
mascons  were  aligned  with  the  sensing  line-of-sight. 

This  test  case  attempted  to  determine  whether  the  mascon's  location 
and  the  sensing  geometry  affected  the  estimated  gravity  field.  The  non- 
dimensional  fit  residuals7  rms  contradicted  the  global  radial  acceleration 
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analysis.  The  limb  mascon  fit  model  (7.37)  resulted  in  a  better  fit  than  the  face 
mascon  fit  model  (11.80).  These  statistics  could  be  the  result  of  the  near-side 
face  mascon  causing  larger  line-of-sight  accelerations  for  the  orbiting  satellite 
than  the  limb  mascons  produced.  These  statistics  could  also  result  from  the 
measurement's  ability  to  sense  the  near-side  mascon's  disturbance  or  inability 
to  recognize  any  disturbance  from  the  limb  mascons.  Unfortunately,  this 
analysis  does  not  indicate  which  estimated  field  is  better  at  modeling  the 
radial  accelerations  in  the  local  vicinity  of  the  mascons.  Contour  plots  of 
these  radial  accelerations  will  demonstrate  the  fit  models'  ability  to  locally 
model  these  mascon  disturbances. 


7.2.2  Analysis  of  the  Eighth  Degree  and  Order  Fits 

Based  on  the  global  lunar  radial  accelerations  errors,  the  dual  orbiter 
bent  pipe  sensing  scheme  produced  the  best  estimated  gravity  field.  This  case 
had  an  rms  radial  acceleration  error  at  100  km  altitude  of  263.015  milligals. 
The  estimated  fit  to  the  single  orbiter,  earth-based  sensing  scheme  had  radial 
acceleration  errors  with  an  rms  of  320.073  milligals.  The  disparity  between 
the  two  rms  values  provides  a  strong  indication  of  the  advantage  of  including 
lunar  far-side  observations  in  the  estimation  process.  These  errors  also  give 
an  indication  of  how  much  more  difficult  it  was  to  estimate  a  spherical 
harmonic  expansion  in  the  presence  of  79  point  mass  disturbances  than  it  was 
in  the  mascon  test  case  with  just  two  point  masses. 


73  Single  Orbit  State  Uncertainty  Prediction  Analysis 

For  the  next  analysis,  a  low  inclination,  low  altitude  lunar  satellite 
orbit  was  used  to  analyze  the  position  and  velocity  errors  between  the  "true" 
and  estimated  gravity  fields.  The  lunar  gravity  field  truth  model  developed 
in  Section  6.4  and  the  two  8x8  spherical  harmonic  expansion  fit  models 
derived  from  the  two  different  sensing  schemes  were  used  for  this  analysis. 
A  15°  inclination,  100  km  lunar  altitude  satellite  orbit  was  selected  for  this 
analysis.  The  estimated  8X8  spherical  harmonic  fit  models  and  their 
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coefficient  uncertainties  were  first  used  to  predict  the  state  uncertainties  for 
this  orbit.  This  covariance  analysis  was  performed  with  both  the  gravity 
harmonic  coefficient  covariance  matrix  (A*1)  and  this  matrix  multiplied  by 
the  fit  residuals'  rms  (I).  Fit  residuals/  once  again,  are  the  observation 
residuals  divided  by  the  assumed  measurement  errors.  This  evaluation  orbit 
was  then  numerically  integrated  using  the  "true"  lunar  gravity  field.  The 
state  errors  between  the  "true"  orbits  and  those  predicted  from  the  two  8X8 
estimated  spherical  harmonic  models  provided  a  further  measure  of  the  two 
sensing  schemes'  capabilities  and  limitations.  Furthermore,  the  predicted 
uncertainties  were  compared  to  the  errors  to  understand  the  uncertainty 
prediction  accuracy  with  gravity  field  modeling  errors. 

The  equations  of  motion  and  equations  of  the  partial  derivatives  of 
motion  with  respect  to  the  gravity  harmonic  coefficients  were  numerically 
integrated  for  one  orbit  (117.85  minutes)  for  the  two  estimated  fit  models 
from  16  May  1968  (JD  2,440,001.5  0h  CT)  with  the  satellite  selenographic  initial 
osculating  orbital  elements  listed  in  Table  7.3-1  (1950.0  angles  in  parentheses). 
Based  on  these  initial  conditions,  the  lunar  orbiter  has  emerged  from  behind 
the  far  side  of  the  moon  and  is  in  the  middle  of  a  near-side  pass.  In  the 
numerical  integration  of  the  "truth"  model,  partial  derivatives  were  not 
calculated.  Auxiliary  software  programs  then  used  the  numerical  integration 
output  files  for  the  three  cases  as  well  as  the  covariance  files  for  the  two 
estimated  fit  models  to  perform  the  error  analysis. 

Table  73- 1:  Low  Inclination,  Low  Altitude  Evaluation  Orbit 
Satellite  Initial  Conditions 


ao 

1838  km  (1.22862712  x  10*5  AU) 

eo 

0.01 

io 

15°  (30.7285027175189°) 

Qo 

210°  (26.9029402084826°) 

(00 

180°  (222.6750569007985°) 

Mo 

315° 
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7.3.1  State  Uncertainty  Prediction 

The  state  uncertainty  prediction  based  on  PEP's  gravity  harmonic 
coefficient  covariance  matrix  (A*1)  derived  from  the  assumed  measurement 
errors  was  overly  optimistic  for~both  cases  since  the  fit  residuals  for  the  two 
estimated  runs  were  so  high  (365.954  and  361.099  rms).  This  method 
predicted  root  sum  squared  (rss)  state  uncertainties  less  than  500  meters  in 
position  and  one  half  of  a  meter  per  second  in  velocity  over  the  orbit  for  the 
single  orbiter  estimated  fit  case.  For  the  dual  orbiter  estimated  fit  model  this 
method  predicted  rss  uncertainties  of  less  than  250  meters  in  position  and  one 
quarter  of  a  meter  per  second  in  velocity  over  the  same  orbit.  These 
predictions  are  extremely  optimistic/  especially  when  compared  to  the  orbital 
errors  presented  in  Section  7.3.2. 

The  prediction  based  on  the  estimated  fit's  covariance  matrix  £  from 
Equation  (5.2-2)  provided  a  more  realistic  assessment  of  the  state 
uncertainties.  This  method  predicted  rss  state  uncertainties  as  high  as  4 
kilometers  in  position  and  4  meters  per  second  in  velocity  for  the  single 
orbiter  estimated  fit  case.  These  predicted  state  uncertainties,  transformed  to 
the  local  vertical/  local  horizontal  coordinates/  are  plotted  over  time  for  this 
lunar  orbit  in  Figure  7.3.1-1.  For  the  dual  body  estimated  fit  model/  this 
analysis  predicted  rss  uncertainties  of  approximately  750  meters  in  position 
and  0.6  meters  per  second  in  velocity.  These  local  vertical/  local  horizontal 
predicted  state  uncertainties  are  plotted  over  time  in  Figure  7.3.1-2.  These 
state  uncertainty  plots  show  that  the  largest  position  uncertainties  are  in  the 
range  direction  and  the  largest  velocity  uncertainties  are  in  the  vertical 
direction  as  might  be  expected  for  the  radial  acceleration  perturbations  caused 
by  the  ma scons. 

This  state  uncertainty  prediction  analysis  demonstrates  the  importance 
of  including  lunar  far-side  observations  in  the  estimation  process.  The 
uncertainties  predicted  for  the  gravity  model  estimated  from  the  dual  orbiter 
sensing  scheme  were  approximately  one  fourth  of  those  estimated  from  the 
single  orbiter  sensing  scheme.  The  dual  orbiter  cases'  uncertainties  were  also 
very  low  while  the  lunar  satellite  passed  across  the  lunar  near  side  (first  forty 
minutes  of  numerical  integration).  For  both  cases,  the  uncertainties  grow  as 
the  orbiter  passes  behind  the  far  side  of  the  moon,  reaching  local  maximums 
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for  velocity  in  the  vertical  direction  and  position  in  the  down  range  direction 
in  the  middle  of  the  far-side  pass.  These  uncertainties  subside  when  the 
spacecraft  reemerges  from  the  far  side  of  the  moon  eighty-eight  minutes  into 
the  orbit.  From  this  point  the  uncertainties  increase  during  the  final  lunar 
near-side  pass. 


7.3.2  State  Errors  Between  "True*  and  Estimated  Gravity  Models 

The  errors  between  the  "true"  and  estimated  fit  lunar  gravity  models 
were  determined  from  a  direct  comparison  of  the  two  numerical  integration 
files.  These  state  errors  were  then  transformed  to  local  vertical,  local 
horizontal  coordinates.  The  local  vertical,  local  horizontal  state  errors 
between  the  orbits  generated  for  the  single  orbiter  estimated  lunar  gravity 
field  and  the  "true"  gravity  field  are  plotted  in  Figure  7.3.2-1.  The  same  errors 
resulting  from  the  difference  between  the  dual  orbiter  estimated  gravity  field 
and  the  "true"  gravity  field  are  plotted  versus  time  in  Figure  73.2-2. 

The  errors  for  the  single  orbiter  estimated  gravity  field  stay  relatively 
flat  as  the  orbiter  crosses  the  near  side  of  the  moon.  Forty  minutes  into  the 
numerical  integration  the  orbiter  passes  behind  the  moon.  Midway  through 
the  far-side  pass  there  is  a  rapid  growth  of  velocity  errors.  Shortly  thereafter 
these  velocity  errors  manifest  themselves  as  position  errors  as  large  as  15 
kilometers  in  the  range  and  10  kilometers  in  the  cross  track  directions.  These 
errors  demonstrate  this  estimated  fit  model's  inability  to  estimate  the  far-side 
gravity  field. 

The  errors  for  the  dual  orbiter  estimated  gravity  field  are  much  lower 
than  in  the  previous  case.  When  plotted  on  the  same  scales  as  the  errors  in 
the  first  case,  the  dual  orbiter  case's  velocity  errors  stay  in  a  narrow  band 
around  zero  meters  per  second.  The  dual  orbiter  position  errors  are  also 
much  more  reasonable  for  the  cross  track  and  vertical  directions.  The  dual 
orbiter  estimated  gravity  model  still  results  in  large  range  errors  over  the 
single  lunar  orbit.  Since  the  dual  orbiter  estimated  lunar  gravity  model  is 
based  on  observations  over  the  entire  lunar  surface,  the  position  and  velocity 
errors  do  not  manifest  themselves  at  a  specific  point  in  the  lunar  orbit  as  was 
the  previous  case.  In  this  case  the  cross  track  and  vertical  position  errors  and 
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Figure  73.1-1:  Single  Orbiter  Estimated  Gravity  Field  State  Uncertainties 
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Single  Orbit  Spacecraft  State  Uncertainty 

Predicted  State  Uncertainty  for  Dual  Orbiter  Sensing  Fit 
Inclination  =  15  deg,  Eccentricity  =  0.01 ,  Altitude  *  100km 


Figure  7.3.1-2:  Dual  Orbiter  Estimated  Gravity  Field  State  Uncertainties 
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Single  Orbit  Spacecraft  State  Errors 

Errors  Between  Single  Orbiter  Sensing  Scheme  Fit  and  Truth  Model 
Inclination  =  15  deg,  Eccentricity  =  0.01 ,  Altitude  =  100km 
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Figure  732-1:  Single  Orbiter  Estimated  Gravity  Field  State  Errors 
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Single  Orbit  Spacecraft  State  Errors 


Errors  Between  Dual  Orbiter  Sensing  Scheme  Fit  and  Truth  Model 
Inclination  =  15  deg,  Eccentricity  =  0.01 ,  Altitude  =  100km 


Figure  7. 3.2-2:  Dual  Orbiter  Estimated  Gravity  Field  State  Errors 
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the  range  and  cross  track  velocity  errors  oscillate  about  zero  for  the  single 
lunar  orbit.  The  vertical  velocity  errors  and  resulting  down  range  position 
errors,  however,  grow  without  oscillating  about  zero  for  the  single  orbit 
propagation. 

For  both  cases  the  state  uncertainty  prediction  was  very  optimistic, 
even  after  including  the  rms  of  the  fit  residuals.  This  discrepancy  results 
from  the  difference  between  the  structure  of  the  two  lunar  gravity  models 
("truth"  and  fit).  The  rss  of  the  predicted  position  uncertainty  was  an  order  of 
magnitude  less  than  the  "true"  position  errors.  The  predicted  uncertainty  in 
velocity  was  only  slightly  better.  Additionally,  the  covariance  analysis  does 
not  successfully  predict  cross  track  position  uncertainties.  This  is  especially 
apparent  for  the  single  orbiter  estimated  gravity  field  case.  Figure  7.3.2-1 
shows  significant  cross  track  state  errors  and  no  notable  cross  track 
uncertainties  are  predicted  in  Figure  7.3.1-1. 

When  analyzing  the  proposed  Lunar  Observer  mission  with  different 
truth  and  fit  models,  Alex  Konopliv  also  noted  that  a  covariance  uncertainty 
analysis  was  unbelievable  because  it  was  "overly  optimistic  for  all  cases"  [25]. 
This  deficiency  in  the  covariance  analysis  is  an  indication  that  the  system 
dynamic  model  needs  to  include  some  process  noise  to  account  for  the  gravity 
field  mismodeling.  A  Kalman  filtering  or  maximum  likelihood  system 
identification  technique  could  include  this  process  noise  as  it  propagated  the 
satellite  equations  of  motion. 


7.4  Lunar  Deorbit  Maneuver  Error  Analysis 

Finally,  a  satellite  lunar  landing  from  a  low  inclination  circular  orbit 
was  simulated  for  both  estimated  gravity  fields.  The  deorbit  bum  for  this 
lunar  maneuver  was  determined  based  on  numerical  orbit  integrations  using 
the  estimated  lunar  gravity  fields.  This  maneuver  was  then  executed  in  the 
"true"  lunar  gravity  field,  and  the  position  errors  at  PDI  were  used  to  evaluate 
the  two  estimated  gravity  field  models'  ability  to  plan  future  lunar  missions. 
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This  analysis  assumed  that  the  spacecraft  was  resupplying  a  lunar  base 
on  the  near  side  of  the  moon  and  had  been  inserted  into  a  circular,  low 
inclination  orbit.  The  transfer  from  this  orbit  to  the  lunar  surface  was 
planned  using  a  Hohmann  transfer.  The  circular  parking  orbit  was 
numerically  integrated  from  the  near  side  of  the  moon  to  the  middle  of  the 
far  side  of  the  moon.  At  this  point,  a  deorbit  bum  placed  the  spacecraft  in  an 
elliptical  transfer  orbit.  The  spacecraft's  pre-planned  Powered  Descent 
Initiation  (PDI)  location  coincided  with  the  elliptical  transfer  orbits'  perilune. 
From  this  point,  the  spacecraft  would  begin  its  powered  descent  to  the  lunar 
surface  and  rendezvous  with  the  lunar  base.  The  powered  descent  portion  of 
the  landing  mission  was  not  simulated  in  this  analysis. 


Figure  74*1:  Lunar  Deorbit  Mission 


The  mission  simulation  began  by  integrating  the  spacecraft's  circular 
orbit  starting  on  16  May  1968  (JD  2,440,001.5  0*0  with  the  satellite  selenographic 
initial  conditions  listed  in  Table  7.4-1  (1950.0  angles  in  parentheses).  Based  on 
the  satellite's  state  at  M=180°  in  the  parking  orbit,  an  initial  guess  for  the  Av 
was  calculated  for  a  Keplerian  transfer  orbit  with  the  PDI  perilune  altitude, 
the  lunar  surface  for  this  analysis.  The  elliptical  transfer  orbit  was  then 
numerically  integrated  and  the  guessed  Av  updated  based  upon  the  PDI 
position  error  at  perilune.  These  iterations  were  continued  until  acceptable 
target  PDI  positions  were  obtained.  Since  this  maneuver  was  planned  for  the 
two  different  estimated  lunar  gravity  fields,  two  different  Ay's  were  calculated 
and  subsequently  executed. 
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For  the  single  body  estimated  gravity  field  case,  the  calculated  deorbit 
bum  occurred  over  0.0611°  south  latitude  and  179.9215°  west  longitude 
approximately  63.985  minutes  into  the  numerical  integration.  The  transfer 
Av  had  a  magnitude  of  2.535  x  10*5  AU/day  (PEP'S  units)  with  a  projected  PDI 
at  0.040°  south  latitude  and  1.951°  west  longitude.  This  burn  was  then 
executed  by  the  spacecraft  in  the  "true"  lunar  gravity  field.  At  63.985  minutes 
into  the  numerical  orbit  integration  the  previously  calculated  Av  was 
subtracted  from  the  satellite's  inertial  velocity  to  simulate  the  deorbit  bum. 
The  numerical  integration  then  proceeded  from  this  new  state.  This 
maneuver  resulted  in  a  PDI  at  0.766°  south  latitude  and  9.903°  west  longitude, 
a  7.952°  error  in  longitude  and  0.725°  error  in  latitude! 

Table  7.4-1:  Low  Inclination,  Lunar  Landing  Parking  Orbit 
Satellite  Initial  Conditions 


ao 

1938  km  (1.295472995  x  10*5  AU) 

eo 

0.0005 

io 

5°  (18.2129222440553°) 

Oo 

0°  (349.1664497093493°) 

«*> 

0°  (226.5049507382479°) 

Mo 

0° 

The  deorbit  was  then  planned  with  the  dual  orbiter  estimated  lunar 
gravity  field.  The  deorbit  bum  was  again  scheduled  for  63.985  minutes  into 
the  orbit  (from  the  satellite  initial  conditions  in  Table  7.4-1)  at  a  selenographic 
latitude  of  -0.0358°  and  longitude  of  -179.9851°.  The  planned  magnitude  of 
the  Av  was  2.5225  x  10-5  AU/day  with  a  projected  PDI  of  0.4339°  south  latitude 
and  5.8691°  west  longitude.  This  maneuver  was  then  executed  with  the 
"true"  lunar  gravity  field.  This  simulation  resulted  in  a  PDI  at  0.2503°  south 
latitude  and  4.008°  west  longitude.  The  PDI  errors  in  this  case  were  now 
1.8613°  in  longitude  and  0.1835°  in  latitude,  a  significant  reduction  in  the  PDI 
position  errors  and  further  evidence  of  the  importance  of  including  lunar  far- 
side  observations  in  lunar  gravity  field  estimations. 

On  the  lunar  surface  these  errors  in  latitude  and  longitude  would 
result  in  a  position  error  of  242  kilometers  in  the  single  orbiter  estimated  8X8 
field  case  and  56  kilometers  for  the  dual  orbiter  estimated  8X8  field.  If  the 
landing  spacecraft  began  powered  descent  this  far  off  target  it  should  still  be 
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able  to  land  dose  to  the  lunar  base,  as  long  as  it  is  dose  enough  to  receive  a 
radio  beacon  signal.  The  suboptimal  descent  trajectory  required  to 
compensate  for  these  errors,  espedally  in  the  first  case,  would  seriously  cut 
into  any  propellant  margins,  perhaps  jeopardizing  future  mission  plans. 

The  magnitude  of  the  navigation  errors  might  be  reduced  if  higher 
than  8X8  spherical  harmonic  or  other  gravity  model  fits  were  estimated.  For 
higher  degree  and  order  fits,  it  would  be  interesting  to  note  whether  the 
above  factor  of  4.3  improvement  in  navigation  accuracy  still  held  when  far- 
side  satellite-to-satellite  measurements  were  added  to  near-side  earth-based 
measurements. 


7.5  Gravity  Coeffident  Parameter  Correlation  Analysis 

Each  of  the  estimation  runs  in  Chapter  Six  produced  extremely  high 
correlations  for  some  of  the  gravitational  coefficients.  These  high 
correlations  inhibited  the  estimation  routine's  ability  to  converge  quickly  and 
provide  an  accurate  estimate  of  the  gravitational  parameters.  Different 
measurement  types  and  orbital  orientations  were  analyzed  to  see  if  any  of 
these  parameter  correlations  could  be  broken.  For  this  analysis,  "truth" 
model  observations  were  generated  for  the  new  observation  methods.  Based 
on  each  new  set  of  observations,  the  lunar  gravity  field  was  estimated  for  a 
single  iteration  of  the  process  outlined  in  Figure  6.1-1.  The  gravitational 
parameter  correlations  were  then  obtained  from  the  single  iteration's 
covariance  matrix. 

For  the  first  correlation  analysis,  the  elliptical  "viewing"  orbit  in  the 
dual  orbiter,  bent  pipe  observation  scheme  was  placed  in  an  orthogonal, 
rather  than  coplanar,  orbit  plane.  The  "sensing"  orbit  was  left  in  its  polar 
orbit  to  provide  full  lunar  surface  coverage  over  the  fourteen  day  observation 
period,  so  the  elliptical  satellite  was  placed  in  an  equatorial  orbit.  From  this 
orbital  geometry,  bent  pipe  observations  identical  to  the  coplanar  case  were 
simulated  for  the  estimation  process.  From  this  iteration,  sixteen  of  the  3,926 
correlations  were  greater  in  magnitude  than  0.90.  Fourteen  of  these  were 
greater  than  0.95.  The  C55C75,  S55S75  pairs  had  the  highest  correlation  of 
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-0.999.  The  C54C-4,  $54574  correlation  pairs  were  the  second  highest  at  -0.9%. 
The  €66^86  and  5 66^86  pairs  also  had  correlations  greater  than  -0.993 

The  second  method  investigated  an  interferometric  observation 
method  to  determine  its  impact  on  the  high  gravitational  parameter 
correlations.  For  this  scenario,  NASA  DSN  stations  made  milli-arc  second 
long  baseline  interferometer  measurements  of  the  near-circular  200  km 
altitude  lunar  orbiter  in  addition  to  the  3  meter  range  and  1  mm/sec  Doppler 
observations.  Because  of  difficulties  encountered  simulating  PEPs  internal 
interferometry  observations,  these  interferometer  measurements  were 
simulated  by  azimuth  and  elevation  angle  observations  for  a  single  DSN  site 
(Goldstone)  with  the  interferometric  accuracy.  Since  the  interferometer 
observations  only  provide  a  planar  angular  measurement,  the  fitting  process 
was  executed  ignoring  the  elevation  angle  observations.  The  fitting  process 
was  then  repeated  using  both  of  the  angle  measurements  to  see  what  impact 
precise  earth-based  three-dimensional  angular  observations  combined  with 
range  and  Doppler  observations  would  have  on  the  gravitational  parameter 
coefficients. 

For  both  of  these  fit  cases  (azimuth  only/azimuth  and  elevation), 
twenty  two  of  the  3,403  correlations  were  greater  in  magnitude  than  0.90  and 
fifteen  of  these  were  greater  than  0.95.  The  parameter  correlations  were 
almost  identical  for  the  two  cases,  with  the  highest  parameter  correlations 
differing  only  in  the  fourth  or  fifth  decimal  place.  The  S55S75  and  C55C75 
pairs  were  the  most  highly  correlated  (-0.999).  The  C65C85,  S65S85,  C54C74, 
554574  parameter  correlations  were  the  second  highest  group.  For  the 
interferometric  case  the  66, 86  correlations  are  still  very  high  (>  -0.990). 

For  the  final  case,  the  Goddard  Space  Flight  Center's  (GSFC)  proposed 
co-orbital  laser  ranging  and  Doppler  lunar  gravity  field  sensing  scheme  was 
simulated  (Figure  1.4-2).  The  main  satellite  was  placed  in  the  polar  near¬ 
circular  200  km  altitude  orbit  which  has  proven  so  useful  throughout  this 
thesis.  The  co-orbiting  subsatellite  was  placed  in  the  identical  orbit  but  with 
an  initial  mean  anomaly  nine  degrees  ahead  of  the  main  satellite.  Laser 
ranging  and  sixty  second  Doppler  observations  were  processed  for  the  entire 
fourteen  day  orbit  with  1  mm  laser  ranging  and  1  mm/sec  Doppler  accuracies. 
These  accuracies  were  based  on  GSFC  briefed  capabilities  [2].  Additionally  this 
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method  used  two-way  coherent  Doppler  observations  between  the  DSN 
Goldstone  site  and  the  main  satellite  using  the  S-Band  (2.115  GHz)  frequency 
with  1  mm/ sec  accuracy  for  sixty  second  count  intervals. 

For  this  case,  the  parameter  correlations  were  different  from  all  of  the 
previous  iterations.  The  highest  gravity  coefficient  correlations  were  now  the 
§54§74  and  C54C74  pairs  (-0.998  and  -0.996  respectively)  with  the  two  initial 
osculating  orbital  eccentricities  also  highly  correlated  (-0.997).  The  66,  86  and 
55,  75  correlations  in  this  case  were  smaller  in  magnitude  than  0.850  except  for 
the  $55375  correlation  which  was  -0.926. 

Each  parameter  estimation  rim  also  determined  the  spread  of  the 
parameter  correlations.  For  each  different  measurement  type  and  orbital 
geometry  investigated,  this  distribution  was  divided  by  the  total  number  of 
correlations.  Figure  7.5-1  plots  this  normalized  distribution  for  the  parameter 
correlations  greater  than  0.50.  This  graph  shows  how  the  single  orbiter,  earth- 
based  sensing  scheme  is  characterized  by  very  high  parameter  correlations.  If 
a  subsatellite  and  laser  ranging  and  Doppler  instrumentation  are  added  to  this 
scheme,  the  correlations  are  driven  down  significantly  and  lunar  far-side 
motion  is  observed.  If  these  mission  modifications  are  not  feasible,  then 
augmenting  the  single  orbiter  earth-based  sensing  scheme  with  long  baseline 
interferometer  measurements  will  also  reduce  the  parameter  correlations, 
although  several  of  the  highest  correlations  still  remain  in  both  cases. 

According  to  the  graph,  a  better  mission  modification  would  be  to 
include  an  elliptical  "viewing7'  satellite  to  observe  the  near-circular  low 
altitude  polar  "sensing"  satellite.  The  chart  suggests  that  if  the  "viewing" 
satellite  is  placed  in  an  orbital  plane  orthogonal  to  the  "sensing"  satellite's 
orbital  plane,  the  lowest  parameter  correlations  are  achieved.  If  both 
spacecraft,  however,  are  launched  together  and  then  separated  after  Lunar 
Orbit  Insertion  (LOI),  as  was  planned  for  the  Lunar  Observer  mission  [39], 
then  coplanar  "viewing"  and  "sensing"  satellites  will  still  provide  one  of  the 
lowest  sets  of  parameter  correlations. 

These  attempts  to  break  the  parameter  correlations  demonstrate  that 
although  some  of  the  parameter  correlations  can  be  reduced,  using  a  near¬ 
circular  polar  200  km  altitude  satellite  as  the  "sensing"  vehicle  to  estimate  an 
8X8  spherical  harmonic  expansion  to  model  the  lunar  gravity  field 
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consistently  results  in  extremely  high  gravitational  parameter  correlations 
between  the  €55 €75,  5 55 3 75,  366386/  and  € 66^86  gravitational  parameters. 
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Figure  7-5-1:  Parameter  Correlation  Distributions  for  Several  Observation  Methods 

Because  the  highest  parameter  correlations  are  always  for  the  n-2,  n-2 
and  n,  n-2  parameters  and  the  n-3,  n-3  and  n-1,  n-3  parameters,  die  spherical 
harmonic  expansion  geometry  should  be  investigated.  Both  sets  of 
correlations  involve  sectorial  terms  and  a  tesseral  counterpart  with  the  same 
number  of  longitudinal  slices.  The  n=m=6  sectorial  and  n=8,  m=6  tesseral 
terms  zero  lines  are  illustrated  by  the  globes  in  Figure  7.5-2.  The  geometrical 
relationship  of  these  highly  correlated  gravitational  parameters  suggests  that 
they  may  result  from  estimating  the  lunar  gravity  field  solely  from 
observations  of  a  polar  lunar  satellite,  since  the  satellite's  ground  track 
repeatedly  traverses  the  sectorial  slices  of  the  moon. 

This  analysis  suggests  that  the  best  way  to  break  the  high  gravity 
coefficient  correlations  would  be  to  use  multiple  inclination  "sensing" 
satellites.  In  his  lunar  harmonic  gravity  analysis,  Alfred  Ferrari  noted  that 
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"...[W]eU  conditioned  estimates  of  the  gravity  harmonics  are  only  achieved 
when  data  from  many  different  inclinations  are  used,  ..."[19].  Multiple 
inclination  observations  could  be  achieved  by  using  new  polar  satellite  data 
along  with  existing  Apollo-era  near  equatorial  satellite  data.  A  better  scheme 
is  to  place  near-circular,  low  altitude  "sensing"  satellites  in  90°  and  45° 
inclination  orbits  and  observe  their  lunar  far-side  motion  with  a  single 
elliptical  "viewing"  satellite.  This  data  could  also  be  combined  with  existing 
Apollo-era  Doppler  tracking  data  when  estimating  the  lunar  gravity  field. 


Figure  75-2:  Tesseral  8,6  and  Sectorial  6,6  Zero  Line  Patterns 
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Conclusions 


8.1  Summary  of  Results 

This  thesis  investigated  the  ability  of  spherical  harmonic  expansion 
estimates  of  the  lunar  gravity  held  to  predict  low  altitude  lunar  orbits 
globally.  These  lunar  spherical  harmonic  expansions  were  estimated  from 
simulated  observations  of  a  near-circular  polar  lunar  satellite.  Several 
different  observation  geometries  and  measurement  methods  were 
investigated;  two  were  used  to  estimate  die  lunar  gravity  held. 

Since  most  of  the  methods  used  to  derive  the  lunar  gravity  held 
employ  Doppler  observations/  a  test  was  performed  to  determine  the  impact 
sensing  geometry  had  on  gravity  field  estimation.  Since  the  Doppler 
observations  measure  velocity  along  the  line-of-sight,  the  geometrical 
orientation  between  the  Doppler  line-of-sight  and  the  vector  between  the 
mascon  and  orbiting  body  will  affect  the  observation's  ability  to  detect  the 
mascon's  presence.  For  the  first  scenario  in  this  test,  a  pair  of  mascons  were 
placed  on  the  lunar  limbs  and  an  eighth  degree  and  order  spherical  harmonic 
expansion  was  estimated.  For  the  second  scenario,  two  mascons  were  placed 
on  the  lunar  near-  and  far-side  faces  and  the  process  was  repeated. 
Unfortunately,  the  test  case  was  unable  to  conclusively  establish  the  impact  of 
viewing  geometry  upon  the  ability  to  detect  and  model  local  mascon 
disturbances. 

A  lunar  gravity  "truth"  model  was  developed  for  this  thesis  which 
combined  a  fifth  degree  and  order  spherical  harmonic  expansion  and  nine 
major  mass  anomalies  or  mascons.  The  seven  most  significant  lunar  maria 
mass  anomalies,  as  identified  by  Wong  et.  al.  [47],  served  as  the  lunar  near- 
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side  mascons.  The  two  significant  far-side  mascons  were  positioned  to  be 
difficult  to  sense  from  earth-based  observations  and  to  balance  the  lunar 
center  of  mass.  This  truth  model,  since  it  contained  mascon  disturbances,  was 
intended  to  be  a  difficult  gravity  field  to  model  with  a  pure  spherical 
harmonic  expansion. 

Using  the  lunar  gravitational  "truth"  model,  the  near-circular  polar 
lunar  satellite  gravity  field  "sensing"  orbit  was  numerically  integrated  for 
fourteen  days,  half  of  a  lunar  period,  to  provide  complete  lunar  surface 
coverage  as  the  moon  rotated  under  the  orbital  plane.  Since  this  satellite  does 
not  actually  measure  the  lunar  gravity  field,  different  techniques  were  used  to 
observe  this  satellite's  motion  and  estimate  the  gravity  field.  Estimates  of  the 
coefficients  in  spherical  harmonic  expansions  sought  to  match  the 
observations  of  the  "sensing"  satellite's  motion.  Based  on  the  observation 
residuals  between  the  initial  parameter  guesses'  theoretical  observations  and 
the  "truth"  model  observations,  the  initial  conditions  and  harmonic 
coefficients  were  adjusted  and  the  process  was  iterated  until  the  best  fit  to  the 
observations  was  obtained.  Although  different  degree  and  order  spherical 
harmonic  expansion  fits  were  attempted,  the  iterative  process  failed  to 
converge  for  those  above  eighth  degree  and  order,  probably  because  of  the 
difference  between  the  "truth"  (spherical  harmonics  plus  mascons)  and  fit 
(spherical  harmonic  expansion  only)  models. 

Estimated  lunar  gravity  fields  were  obtained  for  two  different  sensing 
schemes.  The  first  scheme  recreated  the  gravitatonal  sensing  method  used 
during  the  Apollo  era  with  mostly  near-equatorial  satellites.  Earth-based 
Doppler  observations  of  the  near-circular  polar  lunar  satellite  were  used  to 
estimate  the  lunar  gravity  field.  The  second  sensing  scheme  employed  a 
second  lunar  satellite  in  an  elliptical  orbit,  viewing  the  first  satellite's  motion. 
Earth-  and  satellite-based  Doppler  observables  simulated  the  coherent  bent 
pipe  link  between  an  earth  tracking  station,  the  elliptical  "viewing"  satellite, 
and  the  circular  "sensing"  satellite  proposed  for  the  Lunar  Observer  mission 
[16,39,40j. 

For  the  cases  in  which  the  estimation  process  converged  upon  a  fit  to 
the  observations,  the  fit  model  was  analyzed  to  determine  how  well  it 
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modeled  the  "true"  lunar  gravity  field.  The  model's  fit  to  the  observations  as 
well  as  its  ability  to  model  global  radial  accelerations  were  investigated. 

The  observation  techniques  which  employed  Doppler  observations  of 
the  lunar  far-side  motion  provided  a  much  better  fit  to  the  "true"  gravity 
field.  Both  estimation  techniques  employed  an  eighth  degree  and  order 
spherical  harmonic  expansion.  In  neither  case,  however,  did  the  eighth 
degree  and  order  spherical  harmonic  expansion  closely  model  the  "true" 
lunar  gravity  field. 

Parameter  covariance  information  determined  in  the  estimation  fitting 
process  was  used  to  predict  satellite  state  uncertainties  ahead  one  lunar  orbit. 
A  low  altitude,  low  lunar  inclination  orbit  was  selected  for  this  analysis. 
Significantly  lower  state  covariances  were  predicted  for  the  estimated  lunar 
gravity  field  based  on  the  dual  orbiter  sensing  scheme.  Observations  of  the 
lunar  polar  satellite  during  far-side  passes  significantly  reduced  the  predicted 
uncertainties  for  the  estimated  fit  model.  In  the  best  case,  the  eighth  order 
spherical  harmonic  expansion  predicted  uncertainties  of  close  to  three 
quarters  of  a  kilometer  in  position  and  one  meter  per  second  in  velocity  for 
the  orbit,  hardly  acceptable  for  future  lunar  missions. 

Using  the  lunar  gravity  "truth"  model,  the  satellite's  "true"  position 
and  velocity  were  numerically  integrated  ahead  one  orbit.  The  true  spacecraft 
state  was  then  compared  to  the  predicted  state  to  reveal  die  state  errors.  These 
errors  were  then  compared  to  the  covariance  uncertainties.  Once  again,  the 
estimated  fit  model  which  included  lunar  far-side  observations  did  a 
significantly  better  job  of  predicting  the  correct  satellite  state,  although  it  still 
provided  errors  as  large  as  2.8  kilometers  in  position  and  2.3  meters  per 
second  in  velocity.  In  neither  case  did  the  spherical  harmonic  expansion  fit 
model's  covariance  uncertainties  come  close  to  predicting  the  correct  state 
errors,  further  evidence  of  the  eighth  degree  and  order  expansion's  inability 
to  successfully  predict  lunar  orbits  for  future  lunar  missions. 

The  estimated  fit  models  were  then  used  to  plan  a  lunar  landing 
deorbit  maneuver.  Starting  on  the  near  side  of  the  moon  in  a  circular  low 
lunar  inclination  parking  orbit,  the  satellite's  orbit  was  numerically  integrated 
to  the  far  side  of  the  moon  where  a  deorbit  maneuver  placed  it  in  an  elliptical 
transfer  orbit.  When  the  spacecraft  reached  a  specified  lunar  altitude,  as 
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detected  by  an  on-board  radar  altimeter,  it  began  its  powered  descent  to  the 
lunar  surface.  After  planning  the  deorbit  maneuver  to  reach  a  specific 
location  for  powered  descent  initiation  (PDI),  the  deorbit  maneuver  was 
executed  in  the  "true"  lunar  gravity  field. 

The  lunar  deorbit  maneuver  planned  with  the  single  orbiter  earth- 
based  Doppler  sensing  scheme  estimated  gravity  field  reached  PDI  altitude 
eight  degrees  in  selenographic  longitude  and  three  quarters  of  a  degree  in 
latitude  away  from  its  target  position,  a  surface  error  of  over  two  hundred  and 
forty  kilometers.  The  lunar  deorbit  maneuver  planned  for  the  dual  orbiter 
gravity  field  was  only  one  and  three  quarters  of  a  degree  of  longitude  and  one 
fifth  of  a  degree  of  latitude  off  target  when  it  reached  PDI  altitude.  This 
scenario's  lunar  surface  error  was  fifty-six  kilometers,  an  error  much  easier  to 
recover  from  during  powered  descent  Both  cases,  however,  indicate  that 
further  navigation  aides  will  be  required  during  lunar  operations,  such  as 
radio  beacons  at  a  lunar  base  or  a  global  lunar  navigation  system. 

Extremely  high  parameter  correlations  were  encountered  when 
estimating  the  gravity  coefficients.  Three  different  measurement  schemes 
were  investigated  to  see  if  they  could  reduce  these  parameter  correlations.  In 
the  first  case,  the  elliptical  "viewing"  satellite  was  placed  in  an  orbital  plane 
orthogonal  to  the  circular  "sensing"  satellite's  orbital  plane.  In  the  second 
case,  earth-based  interferometer  measurements  were  added  to  the  single 
orbiter,  earth-based  Doppler  sensing  scheme.  Both  east-west  and  north-south 
look  angles  were  simulated  in  this  search  for  a  way  to  break  the  high 
parameter  correlations.  The  third  case  simulated  the  Goddard  Space  Flight 
Center's  proposed  co-orbital  laser  ranging  and  Doppler  sensing  scheme  [2]. 

The  first  method  was  able  to  reduce  some  of  the  high  parameter 
correlations  from  the  dual  orbiter  co-planar  sensing  scheme,  but  several  of 
the  highest  parameter  correlations  remained.  The  second  case  significantly 
reduced  the  high  parameter  correlations  from  the  single  satellite  earth-based 
Doppler-only  sensing  scheme,  although  once  again  the  same  highest 
parameter  correlations  still  remained.  The  third  case  produced  some 
different  high  parameter  correlations.  This  method  reduced  the  parameter 
correlations  significantly  from  the  single  orbiter,  earth-based  Doppler  sensing 
scheme.  Although  this  case  would  reduce  the  parameter  correlations  from 
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the  single  orbiter  case,  it  does  not  reduce  the  correlations  as  much  as  the  dual 
orbiter  coplanar  and  orthogonal  cases. 

Since  none  of  these  new  observation  techniques  eliminated  the  highest 
parameter  correlations,  the  spherical  harmonic  expansion  geometry  was 
investigated  to  determine  why  the  same  parameters  were  so  highly 
correlated.  It  seems  that  the  extremely  high  correlations  result  from  using  a 
single  polar  "sensing"  orbit  for  all  of  the  gravity  field  estimates.  Polar  ground 
tracks  made  it  difficult  to  separate  the  effects  of  some  of  the  sectorial  and 
related  tesseral  terms.  "Sensing"  orbits  at  different  orbital  inclinations  should 
break  the  extremely  high  parameter  correlations  and  aid  the  lunar  gravity 
field  estimation  process. 


8.2  Recommendations  for  Future  Research 

The  research  for  this  thesis  revealed  several  topics  for  further 
investigation.  First,  a  more  efficient  gravity  field  model  based  on  a  surface 
layer  representation  should  be  investigated.  A  surface  layer  representation 
would  require  roughly  one  third  of  the  parameters  than  a  spherical  harmonic 
expansion  to  model  the  lunar  gravity  field's  high  frequency  behavior.  The 
surface  layer  representation  should  be  able  to  constrain  the  total  lunar  mass 
and  lunar  center  of  mass.  Estimated  fits  with  this  gravity  field  model,  if  it  is 
programmed  into  the  Planetary  Ephemeris  Program,  can  be  compared  to 
those  obtained  with  a  spherical  harmonic  expansion  gravity  field  model. 

Additionally,  this  thesis  only  compared  the  estimated  fit  models  for 
two  observation  techniques.  Because  of  its  observability  into  the  lunar  far- 
side  motion,  the  dual  orbiter  scheme  naturally  provided  a  better  fit  to  the 
lunar  gravity  "truth"  model.  Estimated  fit  models  should  be  developed  for 
the  three  different  observation  methods  investigated  to  break  the  high 
parameter  correlations.  Different  Doppler  observation  methods,  besides  the 
bent  pipe  method,  can  be  investigated,  including  the  satellite  bounce  and 
satellite  beacon  methods  proposed  by  JPL  [40].  Straightforward  software 
changes  in  PEP  should  allow  simulation  of  these  observables,  and  further 
modifications  would  allow  the  incorporation  of  lunar  navigation  aids  (either 
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lunar  surface  sites  or  navigation  satellites).  Additionally,  two  "sensing" 
satellites  in  different  inclinations  should  be  studied  to  verify  their  ability  to 
break  the  high  parameter  correlations.  Higher  degree  and  order  simulated  fits 
should  also  be  attempted  to  determine  the  navigation  improvements  that  can 
be  achieved  with  higher  resolution  spherical  harmonic  expansions. 

For  further  studies,  a  less  stringent  lunar  gravity  "truth"  model  should 
be  used.  The  difficult  to  observe  mascon  on  the  lunar  far  side  was  too  large, 
in  both  surface  area  and  total  strength.  The  lunar  center  of  mass  should  still 
be  constrained,  but  the  far-side  mascons  should  be  roughly  the  same  size  and 
strength  as  their  near-side  counterparts. 

If  lunar  rotation  moment  of  inertia  partial  derivatives  in  the  Planetary 
Ephemeris  Program  are  changed  to  equivalent  second  harmonic  partial 
derivatives,  PEP  can  include  lunar  laser  corner  reflector  observations 
simultaneously  with  satellite  observations  in  the  estimation  of  the  lunar 
gravity  field.  The  lunar  laser  observations  would  provide  very  accurate 
determinations  of  the  lower  degree  harmonics  (second  and  even  third  and 
fourth  degree  [13]),  while  the  lunar  satellite  observations  would  aid  the 
estimation  of  higher  degree  harmonics. 

Finally,  the  most  appropriate  method  to  estimate  the  lunar  gravity 
field  with  real  observations  is  to  use  maximum  likelihood  system 
identification.  This  method  runs  a  Kalman  filter  on  the  satellite  motion  with 
noise  in  the  dynamics  due  to  unmodeled  forces,  and  applies  a  maximum 
likelihood  estimator  to  gravity,  initial  condition,  and  other  parameters  [31 
Chapter  10]. 
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Evaluation  of  Legendre  Polynomials 

and 

Normalized  Legendre  Functions 


A.1  Legendre  Polynomials 

The  Legendre  polynomials  and  their  first  and  second  derivatives  are 
required  for  the  numerical  integration  of  the  equations  for  satellite  motion 
and  for  the  partial  derivatives  of  satellite  motion.  The  Legendre  polynomials 
are  used  to  determine  the  zonal  harmonic  gravity  effects  and  are  also  used  in 
the  recursive  formulas  for  the  Legendre  functions.  Unnormalized  zonal 
harmonic  gravity  coefficients  are  input  to  PEP  since  its  algorithms  use  the 
unnormalized  form  of  the  Legendre  polynomials.  These  algorithms  were  left 
unchanged,  because  the  normalization  factor  in  the  equation 

P„  (z)  =  V2n  +  1  P„  (z)  (A.l-1) 

does  not  grow  too  rapidly  with  zonal  harmonic  degree  n,  and  the  normal 
equation's  automatic  scaling  feature  works  around  any  numerical  problems 
(Section  5.1.6).  A  switch  to  normalized  zonal  harmonic  coefficients  Cn0  from 
Jn  would  also  involve  a  sign  change  in  several  PEP  subroutines. 

The  recursive  evaluation  of  Pn(z),  Pn(z),  P„(z)  in  subroutine  LEGNDR 
(and  the  new  subroutine  LEGNDS)  uses  die  relations  [1, 8] 

nP„  (z)  =  (2  n  - 1  )z  P„_,  (z)  -  (n  -  1)P„_2  (z) 

P'  (z)  =  P'_2  (z)  +  (2n  -  1)P„_,  (z) 

P? (z)  =  P'U  (z)  +  (2n  -  DP'.!  (z) 
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with  the  starting  values 

P0(z)  =  l,P1(z)  =  z  ,P2(z)  =  %(3z2  -1)  ,P3(z)  =  y2(5z3  -3z)  (A.l-5) 

P'0(z)  =  0,  P[{ z)  =  l,  P'2( z)  =  3z,  P3/(z)  =  K(15z2-3)  (A.l-6) 

Pq(z)  =  0 ,  P[\z)  =  0,  P2 (z)  =  3 ,  P3"(z)  =  15z  (A.l-7) 

The  evaluation  of  the  above  functions  starts  with  n  =  2,  since  the 
spherical  harmonic  expansions  in  PEP  start  at  the  second  degree  (assuming 
the  center  of  mass  of  the  central  body  coincides  with  the  origin  of 
coordinates). 


A.2  Normalized  Legendre  Functions 


The  Legendre  functions  and  their  first  and  second  derivatives  are 
required  for  the  numerical  integration  of  the  equations  for  satellite  motion 
and  for  the  partial  derivatives  of  satellite  motion.  Legendre  functions  are 
used  tc  determine  the  tesseral  harmonic  gravity  effects.  The  original  version 
of  PEP  converted  normalized  tesseral  harmonic  gravity  coefficients  (£nm, 
5nm)  from  the  input  stream  into  unnormalized  coefficients  (Cnm,  Snm)  since 
its  internal  algorithms  computed  the  unnormalized  versions  of  the  Legendre 
functions.  The  option  was  added  to  PEP  (incorporated  into  the  SAO  version) 
to  use  normalized  Legendre  functions  in  the  numerical  integration  process. 
This  option  was  desirable  since  the  normalization  factors  vary  widely  for  high 
degrees  n,  especially  as  m  approaches  n  in  the  equation 


2(2rt+l)(n-m)l 


(n  +  m)  I 


P__(z),  m  =  l,...,n  (A.2-1) 


Define  the  precalculated  coefficients  for  m  =  l,...,n  and  n  £  1 


a 


nm 


( 2n  +  l)(n-m ) 
(2n-l)(n  +  m) 


(A.2-2) 
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Km  = 


(2w)(2w  +  l) 
'  (n  +  1) 


m  =  1 


,  (2n  +  l)(n  +  m-l)  . 

b„m  =  -  — - — - —  ,  1  <  m  <  n 

nm  1  (2n-l)(n  +  m) 


(A.2-3) 


=  V(2n)(2n  +  l)(n  +  l)  ,  m  =  1 

-  (A.2-4) 

dnm=^(n  +  m)(n-m  + 1),  l<m<n 

and  the  following  expression  from  the  recursive  formulas  for  the 
unnormalized  Legendre  functions 


q  =  Vl-z2 


(A.2-5) 


The  recursive  evaluation  of  P'nm(z)  is  performed  in  the  new 

subroutine  LEGNDS  using  the  relations  below  [1,  8].  The  recursive  formulas 
for  the  normalized  Legendre  functions  are 

Fu(z)  =  V3i j,  P21(z)  =  i/lS  flz,  P22(z)  =  KVT592  (A.2-6) 

P«l(2)  =  «nl  Z  F„-u(Z)  +  b.l  9  *Vl(Z)  '  ">2  (A-2-7) 

=  zF-taW+i1-  ?  m  =  2,...,n-l  (A.2-8) 

F»„(z)  =  9  (z)  (A.2-9) 

The  recursive  formulas  for  the  first  derivative  of  the  Legendre  functions  with 
respect  to  the  argument  z  (using  unnormalized  Legendre  polynomials)  are 

P,',  (z)  =  ,  P2',(z)  =  — (l-2z2),  P^(z)  =  -zVl5  (A.2-10) 

9  9 


P«'i(z)  =  £(z  P„,  (z)  -  d„,  qPn  (z)) 


(A.2-11) 


pL^  =  ^{mzPm,(z)-dmiqPnim.l(z)),  m  =  2,...,n  (A.2-12) 

Lastly,  the  recursive  formulas  for  the  second  derivative  of  the  Legendre 
functions  with  respect  to  the  argument  z  (using  unnormalized  Legendre 
polynomials)  are 
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4.  P?,(z)  =  ?M{2z2-3),  P£(z> 


=  -Vl5 


p;.(z)=4- 

H 


ii^ip^zHzP^z) 

•f 

-d«{ip-iz)  +  1p’.w) 


(A.2-13) 


(A.2-14) 


f’i(2)  =  4 


^■(l  +  z2)P„„(z)  +  mzP;m(2) 

-C  (f  (2)  +  9  ?»,»-!  (z)) 


,  m  =  2,...,n 

(A.2-15) 


The  recursive  evaluation  of  the  above  functions  starts  with  n  =  2  for 
the  gravity  harmonic  force  evaluation,  since  the  spherical  harmonic 
expansions  in  PGP  start  at  the  second  degree  (assuming  the  center  of  mass  of 
the  central  body  coincides  with  the  origin  of  coordinates). 
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List  of  Constants 


Primary 

Constant 

Description 

Value 

Units 

[201 

M, 

Mass  of  the  sun 

1 

•Jgm, 

Gaussian  gravitational  constant 

0.01720209895 

AIJ3/2/ day 

CT 

Coordinate  Time 

A.1  +  32.15s 
TAI+32.1848 

day 

A.1,  TAI 

Atomic  Time  Seconds 

9,192,631,770 

oscillations 
of  cesium 

The  A.1  and  TAI  atomic  times  are  kept  as  the  average  of  a  number  of  cesium  atomic 
clocks  at  the  national  time  services  (particularly  at  the  U.S.  Naval  Observatory).  UTC  time 
runs  at  the  A.1  and  TAI  rates,  and  presently  differs  from  TAI  by  an  integral  number  of  seconds. 
As  required  (approximately  once  a  year),  there  is  an  increase  of  one  second  (a  UTC  leap  second) 
in  TAI-UTC,  to  keep  UTC  within  one  second  of  UT1  time,  defined  by  a  formula  in  terms  of 
sidereal  time  [23].  Given  a  UTC  observation  time,  CT  is  computed  for  interpolating  from 
ephemeris  files,  and  UT1  time  and  earth  wobble  are  computed  from  tables  published  by  the 
International  Earth  Rotation  Service  (see  Appendix  C.2). 
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Auxiliary 

Constant 

123] 

Description 

Value 

AULTSC 

Astronomical  Unit 

in  Light  Seconds 

499.004782 

c 

Speed  of  light 

299,792.458 

m5 

Mc 

sun/(earth  +  moon)  mass  ratio 

328,900.1 

Mc 

M| 

(earth  +  moon) /moon  mass  ratio 

82.301 

GM/ 

Lunar  gravitational  constant 
(PEP  uses  AtP/day2) 

4,902.79375 

Pi 

Lunar  mean  equatorial  radius 

1,738 

Pi 

Lunar  period 

27322 

gal 

Unit  of  Acceleration 

1.0 

Units 

seconds/AU 

km/sec 

km3 /sec2 

km 

days 

cm/sec2 
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Appendix  C 

Inertial  Coordinate  Transformations 


Cl  Transformation  Between  Moon  Fixed  and  Inertial  Frames 

The  transformation  between  moon-fixed  (selenographic)  coordinates 
and  the  inertial  integration  frame  coordinates  (the  mean  equinox  and  equator 
of  1950.0)  is  required  to  evalute  gravity  harmonic  and  mascon  accelerations  in 
the  numerical  integration  of  a  lunar  satellite's  motion.  Auxiliary  software 
also  uses  this  transformation  to  analyze  integration  output  in  the 
selenographic  frame  and  to  prepare  numerical  integration  input.  In  Equation 
(3.2-2)  for  the  transformation  between  coordinates  fixed  in  the  moon  and  the 

4 

integration  frame  coordinates,  the  rotation  matrix  R  can  be  expressed  as  [6] 

R  =  UVP  (C.l-1) 

where 

P  -  Earth  precession  matrix  transforming  between  integration 
frame  coordinates  and  coordinates  referred  to  the  mean 
equinox  and  equator  of  date  (50  arcseconds  per  year). 

V  =  Transformation  between  coordinates  referred  to  the  mean 
equinox  and  equator  of  date  and  coordinates  referred  to  the 
mean  equinox  and  ecliptic  of  date  (23.4°  rotation). 

U  =  Transformation  between  coordinates  referred  to  the  mean 
equinox  and  ecliptic  of  date  and  coordinates  fixed  in  the 
moon  along  the  nominal  principal  moment  of  inertia  axes 
(360°  per  27.2  days). 
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PEP  uses  a  Taylor  series  expansion  to  calculate  the  earth  precession 
matrix  (P)  for  speed  during  a  numerical  integration  [7].  If  e0  is  the  standard 
expression  for  the  obliquity  of  the  ecliptic  [20,  23],  the  transformation  matrix 
between  mean  equatorial  and  ecliptic  coordinates  of  date  is 

'10  0 
V=  0  cos£0  sine0 
0  -sin £0  cos£q 

Cassini's  laws  plus  the  physical  libration  of  the  moon  determine  the 
transformation  between  the  mean  equinox  and  ecliptic  of  date  and 
selenographic  coordinates  (U).  The  following  formulas  are  used  within  PEP 
to  calculate  this  transformation  matrix. 

Letting 

M  =  Mean  longitude  of  the  moon  measured  in  the  ecliptic  from 
the  mean  equinox  of  date  to  the  mean  ascending  node  of 
the  lunar  orbit  and  then  along  the  orbit  (27.2°  day  period). 

Cl  =  Longitude  of  the  mean  ascending  node  of  the  lunar  orbit 
on  the  ecliptic  measured  from  the  mean  equinox  of  date 
(18.6  year  period). 

I  =  Inclination  of  the  lunar  equator  to  the  ecliptic  (1.53889°). 

o,  p,  x  =  Physical  librations  in  node,  inclination,  and  longitude. 

Standard  polynomial  expressions  are  used  for  the  angles  M  and  £2,  as  well  as 
for  the  angles  £,  V,  F,  and  D  which  are  from  Brown's  lunar  theory  [23].  PEP 
uses  the  following  trigonometric  series  for  the  physical  librations  [28] 

x  =  - 12.9"  sin  £  -  0.3"  sin  2£  +  65.2"  sin  £'  +  9.7"  sin  (2f  -  2£) 

+  1.4"  sin  (2F  -  2D)  +2.5"  sin  (D  -  t)  -  0.6"  sin  (2D  -  2£  +  £') 

-  7.3"  sin  (2D -20  -  3.0"  sin  (2D  -  £)  -  0.4"  sin  2D 

+  7.6"  sin  Cl  (C.l-3) 

p  =  - 106"  cos  £  +  35"  cos  (2F  -  l )  -  11"  cos  2F 

-  3"  cos  (2F  -  2D)  -  2"  cos  (2D  -  l )  (C.l-4) 
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Appendix  G  Inertial  Coordinate  Transformations 


I(x-o)  =  108" sin  /  - 35" sin (2F -  /)  +  ll"sin2F 

+  3"  sin  (2F  -  2D)  +  2"  sin  (2D  -  /)  (C.l-5) 

where  I  is  measured  in  radians. 

In  terms  of  these  angles,  die  Euler  angles  for  die  moon  rotation  are 
y  =  ft  +  <j 

0  =  I+p  (C.l-6) 

*  =  180°+(M-Q)+(t-<7) 

The  transformation  matrix  between  stenographic  coordinates  fixed  in  the 
moon  and  those  referred  to  the  mean  equinox  and  ecliptic  of  date  is  then  [6] 

Un  =  -  sin  y  cos  0  sin  $  -  cos  y  cos  $ 

U|2  =  cosy  cos0  sin$  +  siny  cosy 

U)3  *  -  sin  0  sin  y 

U21  =  -  sin  y  cos  0  cos  $  -  cosy  sin  y 

U22  =  cos  y  cos  0  cos  $  -  sin  y  sin  $ 

U23  =  -  sin  0  cosy 

U31  =  -  sin  y  sin  0 

U32  =  cosy  sin 0 

U33  =  cos  0  (C.l-7) 

To  check  these  transformations/  the  earth  and  sun's  selenographic 
longitude  and  latitude  were  printed  out  for  certain  dates.  These  positions 
were  obtained  by  transforming  the  vectors  pointing  from  the  moon  to  these 
objects  with  the  transformation  matrix  Equation  (C.l-1).  This  check  verified 
that  PEP's  values  agreed  with  the  published  values  in  the  Astronomical 
Almanac,  to  the  number  of  places  published  (0.001°  for  die  earth,  0.01°  for  the 
sun). 


Since  this  transformation's  moon-fixed  axes  are  not  exacdy  along  the 
moon's  principal  moments  of  inertia  axes,  the  C21,  S21,  and  S22  spherical 
harmonic  coefficients  cannot  be  assumed  to  be  zero.  If  the  rotation  and  orbit 
of  the  moon  were  estimated  by  fitting  to  lunar  laser  comer  reflector  data, 
these  second  degree  harmonic  coefficients  could  be  set  to  zero. 
Simultaneously  processing  lunar  laser  data  with  lunar  orbiter  data  could 
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provide  the  best  estimate  of  the  remaining  second  degree  harmonic 
coefficients,  as  well  as  third  and  perhaps  fourth  degree  coefficients. 

One  auxiliary  program  (selenang),  converts  satellite  osculating  orbital 
angles  in  the  selenographic  coordinate  frame  to  angles  in  the  integration 
(1950.0)  frame.  These  angles  are  input  to  PEP  to  specify  a  desired 
selenographic  orbital  orientation  for  numerical  integrations.  Within  PEP, 
these  initial  osculating  elliptic  orbital  elements  are  converted  to  Cartesian 
position  and  velocity  initial  conditions  for  the  numerical  integrations. 

Another  auxiliary  program  (selenelm)  converts  the  numerical 
integration  output  to  selenographic  position  and  velocity  and  to 
selenographic  osculating  elliptic  orbital  elements  as  functions  of  time.  The 
selenographic  osculating  elliptic  orbital  elements,  the  altitude  above  the  lunar 
surface,  and  the  subsatellite  selenographic  longitude  and  latitude  can  then  be 
plotted  over  time. 


C2  Transformation  Between  Earth  Fixed  and  Inertial  Frames 

Since  earth-fixed  observing  sites  need  to  be  transformed  to  the  inertial 
frame  to  process  observations,  the  transformation  matrix  between  earth-fixed 
coordinates  and  the  inertial  integration  frame  referred  to  the  mean  equinox 
and  equator  of  1950.0  is  also  required.  This  rotation  transformation  matrix  R 
can  be  expressed  as  [8] 

R  =  WSNP  (C.2-1) 


where 

P  =  Earth  precession  matrix  transforming  between  integration 
frame  coordinates  and  coordinates  referred  to  the  mean 
equinox  and  equator  of  date  (50  arcseconds  per  year). 

N  =  Earth  nutation  matrix  transforming  between  coordinates 
referred  to  the  mean  equinox  and  equator  of  date  and 
coordinates  referred  to  the  true  equinox  and  equator  of  date 
(20  arcseconds). 
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Appendix  C;  Inertial  Coordinate  Transformations 


S  -  Earth  sidereal  time  matrix  transforming  between 
coordinates  referred  to  the  true  equinox  and  equator  of  date 
and  coordinates  with  z  axis  along  the  earth  pole  of  rotation 
and  x  axis  in  the  meridian  of  Greenwich  through  the 
rotation  pole  (360°  per  23h  56*"  4.09*). 

W  *  Earth  wobble  matrix  transforming  between  coordinates 
with  z  axis  along  the  earth  pole  of  rotation  and  x  axis  in  the 
meridian  of  Greenwich  and  coordinates  with  z  axis  along 
the  conventional  international  pole  and  x  axis  in  the 
meridian  of  Greenwich  through  the  conventional 
international  pole  (0.3  arcseconds). 

Since  PEP  uses  the  IAU  value  of  the  precession  constant  [23]  in  the 
trigonometric  angles  for  evaluating  the  precession  matrix,  P,  numerical 
integration  results  in  the  1950.0  reference  frame  are  transformations  of 
integration  results  in  the  IAU  J2000.0  reference  frame  [23]. 

Since  observations  in  PEP  are  a  function  of  UTC  time  and  there  is  a 
mathematical  formula  relating  sidereal  time  to  UT1  time  [23],  the 
relationship  between  UTC  and  UT1  time  is  needed  to  determine  the  sidereal 
time  transformation  matrix,  S.  PEP  determines  this  relationship  from  values 
published  by  the  International  Earth  Rotation  Service  (IERS)  based  on  very 
long  baseline  interferometry  and  lunar  and  satellite  laser  ranging 
observations. 


UT1  -  UTC  =  (A.1  -  UTC)  -  (A.1  -  UT1)  (C.2-2) 

The  IERS  also  publishes  wobble  angles  based  on  the  above-mentioned 
observations  which  are  used  to  compute  the  earth  wobble  transformation 
matrix,  W.  Exact  values  of  these  quantities  were  not  required  for  this  thesis' 
simulations. 
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Appendix  D 


Lunar  Gravitational  "Truth"  Model 


Table  D-l:  5X5  Spherical  Harmonic  Coefficient  Portion  of  Truth  Model  (12] 


Harmonic 

xlO^ 

Harmonic 

xlO-6 

J2 

202.431 

C21 

-0.07 

§21 

-0.00 

C22 

34.49 

S22 

0.03 

mmmm 

1  iw  ha  a  1  w  ii'iiiii 

S31 

6.63 

c32 

§32 

4.76 

C33 

15.87 

S33 

-2.45 

mm  ■ 

-11.73 

-432 

s« 

ma 

-8.13 

§42 

I  K 

0.48 

S43 

-14.43 

C44 

-330 

S44 

-035 

J5 

2.388 

C51 

-9.66 

§51 

-133 

C52 

3.71 

S52 

-235 

C53 

-039 

S53 

4.91 

C54 

0.56 

S54 

-638 

C55 

-6.69 

S55 

11.60 
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Table  D-2:  Lunar  Near-Side  Masoons  at  1638  km  Radius  [47] 


Mass 

East 

North 

Mass 

East 

North 

(10"6  Lunar 

Longitude 

Latitude 

(10"6  Lunar 

Longitude 

Latitude 

Mass) 

(degrees) 

(degrees) 

Mass) 

(degrees) 

(degrees) 

|  Sea  of  Serenity 

[  Sea  of  Rains 

10.5 

34.0 

342.0 

48.0 

2.437 

17.5 

34.0 

0.467 

3285 

405 

2.453 

105 

28.0 

5.329 

3375 

40.5 

3.316 

175 

28.0 

5.218 

3465 

405 

3.774 

245 

28.0 

1.252 

3285 

34.0 

0.371 

8.25 

22.5 

1.628 

3355 

34.0 

2.628 

!  13.75 

22.5 

3.586 

3425 

34.0 

3.628 

19.25 

22.5 

1.752 

349.5 

34.0  ! 

1.681 

24.75 

22.5 

-0.421 

3285 

28.0 

0.290 

12.5 

17.5 

2.336 

3355 

28.0 

0.624 

175 

17.5 

2.654 

3425 

28.0 

0.229 

225 

17.5 

c 

w 

>ea  of  Crises  | 

0.328 

■TOM 

52.25 

22.5 

1.934 

1 

1 

57.75 

22.5 

0.531 

525 

17.5 

1.522 

3525 

75 

2.837 

57.75 

175 

|  Sea  of  Moisture 

1.836 

62.5 

17.5 

mmmwm 

318.75 

-225 

1.183 

575 

17.5 

t  | 

32525 

-22.5 

0.688 

625 

12.5 

|  0.164 

3215 

-28.0 

Sea  of  Nectar  | 

j  Smyth's  Sea 

0.251 

27.5 

-125 

0.622 

875 

2.5 

1.394 

325 

-125 

-0.172 

82.5 

-25 

1.199 

375 

-125 

0.873 

87.5 

-25 

275 

-17.5 

0.739 

925 

-2.5 

325 

-175 

0.352 

87.5 

-75 

0.663 

375 

-17.5 

0.691 

92.5 

-75 
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Table  D-3:  Lunar  Far-Side  Mascons  at  1638  km  Radius 


Mass 

(1C*  Lunar 
Mass) 


Mass 

(10*  Lunar 
Mass) 


rve  Mascon 


•ftl 


1775 
180.0 
1825 
1.357  176.75 

2.929  178.75 

2.265  181.25 

0.935  183.75 

1.323  175.0 

1775 


ring  Mascon 


-48.0 

-48.0 

-48.0 


Final  Point  Masses 


4.85x10*  214.93  -49.06 


0.829 

182.5 

1.504 

180.0 

1.254 

1775 

0.615 

183.75 

1.391 

181.25 

1.613 

178.75 

0.825 

176.25 

2.125 

185.0 

4.737 

1825 

7.712 

180.0 

3.713 


1.11385  217.6253 

2-lxlO*11  203.77 


kftl 


0.0 

0.0 

0.0 


-505 

-505 

-455 

-455 
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Appendix  E 


Tables  of  Spherical 
Harmonic  Coefficients 


On  the  following  pages  are  the  tables  containing  the  spherical 
harmonic  coefficients  for  the  models  estimated  in  this  thesis.  Also  included 
are  the  spherical  harmonic  coefficients  from  degree  six  through  ten  from  Alex 
Konopliv's  50  X  50  spherical  harmonic  model  estimated  at  the  Jet  Propulsion 
Laboratory.  These  coefficients  were  used  in  the  fit  of  a  10  X  10  model  to 
observations  generated  by  a  10  X  10  truth  model.  The  tesseral  coefficients  are 
normalized  and  the  zonal  coefficients  are  unnormalized  to  follow  the 
convention  used  in  PEP-D  which  was  modified  from  the  SAO's  version  to 
use  normalized  Legendre  polynomials,  but  not  normalized  Legendre 
functions. 
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Table  E-l:  Limb  Mascon  Fit  Model  Coefficients 


Note:  All  of  the  Cnm  and  Snm  terms  are  Normalized.  The  Jn  terms  are  Unnormalized. 
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Table  E-2:  Face  Mascon  Fit  Model  Coefficients 


Harmonic 


212.9685 

-172629 

6.8783 


10.1197 

35.1938 


20.1172 

0.1584 

19.8874 


-25.7774 

-2.1818 

14.6293 

-43549 


-73438 

24.1943 

-8.1387 

-113970 

-4.6557 


6.4501 

-0.9473 


13589 

-9.6129 

2.9190 

10.2397 

-3.1597 

-03898 

0.1187 


-63178 

-2.6772 

3.3579 

-03695 

-5.7318 

3.6244 

-0.4581 

13694 


-18.2734 

567533 

-33.1132 


-53137 


3.4003 

11.9870 

18.4118 

3.9231 


-03299 

-73913 

12.9983 

-07905 

7.4909 


-33481 


iiru!: 


7.6098 

-93866 

-4.1229 

3.4506 


-33029 

4.0109 

-43617 

-4.9557 

7.1339 

1.1869 

0.0073 


-0.4581 

1.4211 

-4.3328 

4.2547 

2.0067 

-53175 

-0.4756 

0.0024 


Note:  AU  of  the  Cnm  and  Snm  terms  are  Normalized.  The  Jn  terms  are  Unnormalized. 
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Table  E-3:  8X8  Single  Orbiter  Earth-based  Doppler  Fit  Model  Coefficients 


Harmonic 


Harmonic 


181.1456 

-74.6007 

133.1974 

-16.5125 


207759 

98.9974 


52.1951 

157778 

-527148 


167472 

-101.0443 

0.9710 

79.4372 


-80.3231 

-6.9808 

137.9127 

18.6358 

-80.5845 


-41.8646 

765292 

0.9371 

-1405093 

-10.4542 

435262 


445845 

157502 

-63.0660 

104385 

125.9903 

245612 

-85265 


15.1214 

-305354 

-6.4764 

33.7568 

-22.1201 

-68.3081 

-55991 

207612 


595185 

-132.1220 

795217 


487268 

635880 


205043 

-22.4302 

-875210 


125355 

-87.6400 

-53.3659 

625998 


117542 

-67.9849 

118.1211 

61.9093 

-35.9216 


-567582 

-9.1320 

765052 

-1157748 

-465469 

245186 


-26.6661 

72.1269 

47952 

-66.4744 

76.7319 

15036 

-05739 


155307 

284645 

-507189 

427387 

365057 

-38.0350 

3.3389 

-1.0150 
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Table  E-4:  Additional  Coefficients  for  the  10  X  10  Spherical  Harmonic  Expansion  12 7} 


Harmonic 


Harmonic 


-0 
-1.5708 
2.0466 
-13200 
-2.8712 
-23382 
-2.6997 
-2.1130 


I  C«J 

■03704 

S99  1 

-3.4853 

mil  m 

-0.4607 

-0.1513 

0.0074 

0.1069 

0.7679 

ClO/4 

-3.4765 

13562 

C103 

13547 

S103 

-0.1451 

Cl  0,6 

0.0773 

Sio^ 

-3.1297 

C107 

-4.6092 

S107 

-1.6688 

Cl0y8 

-23105 

S103 

-6.1743 

Qo^ 

-3.9610 

Sl0,9 

0.0054 

Qo.io 

0.6457 

Sio.io 

2.4016 

Note:  All  of  the  Qim  and  Snm  terms  are  Normalized.  The  Jn  terms  are  Unnormalized. 
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Table  E-5:  8X8  Dual  Orbiter  Bent  Pipe  Doppler  Fit  Model  Coefficients 


Harmonic 


h 


2223067 

163116 

-35341 

62438 


20.9141 

629024 


233231 

173213 

5.0245 


Harmonic 


h 

J5 


255174 

-26.9111 

9.0723 


-165822 

3.6949 

135951 

7.6961 

-28.7475 


8 

-43649 

6.4123 


62918 


-123986 

-23.7569 

27198 


-85788 

-95117 

-0.4007 

-05922 

-33561 


8.9587 

-45860 

85267 

-14.6763 

-142858 

-22788 


5.7153 

-13007 

6.7590 

-135958 

225105 

-1.1883 

0.9170 


3.4853 

-2.9015 

-7.9046 

•02085 

-2.3673 


kAU... 


-03801 

1.0937 


12082 

-1.1762 

-12396 


Note:  All  of  the  Cnm  and  Snm  terms  are  Normalized.  The  Jn  terms  are  Unnormalized. 
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